dist_m4ri.c

Go to the documentation of this file.

 
#define _GNU_SOURCE
#include <inttypes.h>
#include <strings.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <stdatomic.h>
#include <time.h>
#include <unistd.h>
#include <pthread.h>
#include <sched.h>
#include <math.h>
#include <limits.h>
#include <m4ri/m4ri.h>
 
#include "mmio.h"
#include "uthash.h"
#include "util_hash.h"
#include "util_m4ri.h"
#include "util_io.h"
#include "dist_m4ri.h"
#include "dist_cc.h"
 
/* Mutex protecting M4RI's internal non-thread-safe MMC memory cache */
static pthread_mutex_t m4ri_mem_mutex = PTHREAD_MUTEX_INITIALIZER;
 
static inline mzd_t *safe_mzd_from_csr(mzd_t *dst, const csr_t *p) {
  pthread_mutex_lock(&m4ri_mem_mutex);
  mzd_t *res = mzd_from_csr(dst, p);
  pthread_mutex_unlock(&m4ri_mem_mutex);
  return res;
}
 
static inline mzd_t *safe_mzd_init(rci_t r, rci_t c) {
  pthread_mutex_lock(&m4ri_mem_mutex);
  mzd_t *res = mzd_init(r, c);
  pthread_mutex_unlock(&m4ri_mem_mutex);
  return res;
}
 
static inline void safe_mzd_free(mzd_t *M) {
  if (!M) return;
  pthread_mutex_lock(&m4ri_mem_mutex);
  mzd_free(M);
  pthread_mutex_unlock(&m4ri_mem_mutex);
}
 
static inline mzp_t *safe_mzp_init(rci_t length) {
  pthread_mutex_lock(&m4ri_mem_mutex);
  mzp_t *res = mzp_init(length);
  pthread_mutex_unlock(&m4ri_mem_mutex);
  return res;
}
 
static inline void safe_mzp_free(mzp_t *P) {
  if (!P) return;
  pthread_mutex_lock(&m4ri_mem_mutex);
  mzp_free(P);
  pthread_mutex_unlock(&m4ri_mem_mutex);
}
 
static inline double get_time_sec(void) {
  struct timespec ts;
  clock_gettime(CLOCK_MONOTONIC, &ts);
  return (double)ts.tv_sec + (double)ts.tv_nsec * 1e-9;
}
 
/* The splitmix64 output function: a bijective 64-bit hash */
static inline uint64_t mix64(uint64_t z) {
  z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
  z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
  return z ^ (z >> 31);
}
 
static inline uint64_t splitmix64(uint64_t *state) {
  return mix64(*state += 0x9e3779b97f4a7c15ULL);
}
 
static inline int rand_uniform_thread(int max, uint64_t *state) {
  if (max <= 1) return 0;
  return (int)(splitmix64(state) % (uint64_t)max);
}
 
static inline mzp_t * mzp_rand_thread(mzp_t *q, rci_t length, uint64_t *state) {
  if (q == NULL) return NULL;
  for (int i = 0; i <= (int)length - 2; i++) {
    q->values[i] = i + rand_uniform_thread(length - i, state);
  }
  for (int i = length - 1; i < (int)q->length; i++) {
    q->values[i] = i;
  }
  return q;
}
 
#define CC_SPLIT_MAX_LEVELS 3 /* CC work units: at most 3 fixed levels, B^3 units per start index (cc_round_levels()) */
 
/* The CC column-syndrome hash, for the last two levels of a cluster (see cc_walk_last2()).  The hash of a set S of rows
 * of H is the XOR of the random 64-bit values zr[r], r in S (Zobrist hashing): the hash of the syndrome s + h_c, where
 * h_c is the column c of H, is that of s XOR hc[c].  The nonzero columns are sorted by (hash, column), and each
 * distinct hash value `key` has an entry {key, beg, cnt} in the open-addressing table `tab` (linear probing from the
 * slot key >> tshift) for its run of columns cols[beg..beg+cnt-1].  The bit (key & bmask) of the filter `bits` is set
 * for each key, so that most failed lookups only read one bit. */
typedef struct {
  uint64_t key;
  int beg, cnt; /* cnt = 0: empty slot */
} cc_hash_slot_t;
 
typedef struct {
  uint64_t *zr;         /* random value of each row of H */
  uint64_t *hc;         /* hash of each column of H (0 for a zero column) */
  uint64_t *bits;       /* filter: bit (key & bmask) set for each key */
  uint64_t bmask;
  cc_hash_slot_t *tab;  /* 2^(64-tshift) slots */
  int tshift;
  int tmask;
  int *cols;            /* the nonzero columns sorted by (hash, column) */
  int ncols;            /* their number */
  int nkeys;            /* number of distinct keys */
} cc_hash_t;
 
typedef struct {
  params_t *p;
  const cc_hash_t *cc_hash; /* CC column-syndrome hash for cc_walk() (NULL with smax > 0: start_CC_recurs_mt()) */
  int num_threads;
  int rw_threads; /* workers 0..rw_threads-1 may run RW (dense RW matrices); all workers may run CC */
  double timeout;
  double start_time;
  int dexp;
 
  /* Distance bounds & stop flags (cache-line isolated for read-mostly access) */
  _Alignas(64) atomic_int dmin; /* dmin-1 is max cluster size analyzed without success */
  atomic_int dmax;              /* smallest weight codeword found (0 if none) */
  atomic_int cc_found_weight;   /* smallest weight of a codeword found by CC (0 if none) */
  atomic_bool cc_exact;         /* cc_found_weight is certified exact (found in a round w=dmin) */
  atomic_bool stop_flag;        /* signals all threads to terminate */
  atomic_bool rw_stop_flag;     /* signals RW workers to stop (in method 3) */
  atomic_int stop_reason;       /* why the run ended (STOP_NONE until then), see ctx_stop() */
  atomic_int cw_feed_limit;     /* codeword_feed_limit() (`min_hits` statistic), updated under `cw_mutex` */
 
  /* RW state (cache-line isolated from CC and read-mostly bounds) */
  _Alignas(64) long total_rw_steps;
  atomic_long rw_steps_started;
  atomic_long rw_steps_completed;
  atomic_llong rw_busy_ns;           /* thread time in ns of the timed RW steps, for the method-3 timing model */
  atomic_long rw_timed_steps;        /* number of RW steps timed in `rw_busy_ns` */
  atomic_int rw_rank;                /* rank of H from the RW elimination (n - nu with ksub > 0), 0 if not known */
  atomic_long rw_unif_steps;         /* completed full-matrix RW steps with uniform random permutations (no window) */
  atomic_int rw_pairs_l;             /* the number l of pivot rows of the RW pairs (0: none), see rw_pairs_step() */
  atomic_llong rw_pairs_checked;     /* pairs of non-pivot columns checked (the weight of their sum computed) */
  atomic_llong rw_hits1;             /* RW codewords recorded with weight up to the upper bound: from single columns */
  atomic_llong rw_hits2;             /* and from the pairs */
  const uint64_t *rw_lhash;          /* hashes of the columns of L for the RW pairs (NULL: none), see rw_lhash_init() */
 
  /* CC state for current weight (cache-line isolated) */
  _Alignas(64) atomic_llong cc_next; /* weight, fixed levels, next work unit of the current CC round, see `cc_pack()` */
  csr_t *mHT_cc;
  int max_col_W;
  bool even_weight_only;                 /* 1_n in rowspan(H) and smax=0: all codewords have even weight */
  int cc_B;                              /* CC work units: radix B, the maximum row weight of H */
  int cc_R[CC_SPLIT_MAX_LEVELS + 1];     /* work units per start index with L fixed levels: B^L (0 if too large) */
  double cc_branch;                      /* typical number of branches per level: average row weight of H - 1 */
  atomic_int cc_active_workers;
  atomic_int cc_target_workers;
  atomic_int cc_round_active;
  atomic_llong cc_busy_ns;           /* total CC thread time in ns (all rounds), for the CC timing model */
  atomic_llong cc_nodes;             /* total calls of the CC recursion (smax = 0, all rounds), for DBG_TIMING */
 
  /* Codeword synchronization */
  pthread_mutex_t cw_mutex;
  /* irred=1: the reducible codewords not recorded (the union of a codeword and a disjoint nonzero zero-syndrome
   * vector), found by CC and by RW (pairs of disjoint circuits, ksub rows) */
  atomic_llong cw_red_cc;
  atomic_llong cw_red_rw;
 
  /* Subspace sketching (ksub) state */
  mzd_t *N_global;
  int nu;
  pthread_rwlock_t basis_rwlock;
  atomic_long next_refresh_step;
 
  /* Timing stats: CC work of each completed round in thread-seconds (total CC thread time) */
  double cc_time_per_weight[MAX_W];
 
  /* Coordinator only: the current CC round (for the status and stop messages), and the periodic status */
  int round_open;         /* 1 from the start of a CC round until it is completed */
  int round_w;            /* weight of the current (or last) CC round */
  int round_beg;          /* its first start index */
  int round_end;          /* its last start index */
  int round_R;            /* its number of work units per start index */
  int round_threads;      /* its planned number of CC threads */
  double round_t0;        /* its start time */
  double round_est;       /* its estimated CC work in thread-seconds (0: unknown) */
  const char *cc_idle;    /* method 3: why no CC round runs (NULL: none) */
  int stop_w;             /* CC weight for the stop message (STOP_CC_DONE, STOP_CC_SLOW) */
  double status_next;     /* time since the start of the next status line (DBG_STATUS) */
  double status_dt;       /* interval to the status line after that */
  double status_t;        /* time since the start of the last status line */
  long status_steps;      /* RW steps completed at the last status line */
 
  /* Thread handles */
  pthread_t *threads;
  struct worker_arg *wargs; /* the arguments of the num_threads workers */
  int spawned;              /* number of workers started (main thread only), see ctx_spawn() */
  int spawn_initial;        /* number of workers started at first (fewer in method 2, see CC_SPAWN_INITIAL) */
  bool spawn_failed;        /* a worker thread could not be created: no more are started */
} distfork_ctx_t;
 
/* Reasons why the run ended (the first one recorded, see ctx_stop()), reported with DBG_SUMMARY */
enum {
  STOP_NONE = 0,
  STOP_BOUNDS,       /* the bounds coincide: dmin = dmax */
  STOP_CC_FOUND,     /* CC found a codeword (the exact distance, unless lower weights were not scanned) */
  STOP_RW_CONVERGED, /* method 1: RW `min_hits` convergence */
  STOP_RW_STEPS,     /* method 1: all RW steps done */
  STOP_TIMEOUT,      /* timeout reached */
  STOP_CC_DONE,      /* CC done up to the largest weight needed (wmax), and RW ended */
  STOP_CC_SLOW,      /* method 3: the next CC round is predicted not to finish before the timeout, and RW ended */
  STOP_WMIN,         /* a codeword of weight <= wmin found */
  STOP_MAXC,         /* maxC codewords collected */
  STOP_COLLECTED,    /* the CC rounds w = d..d+dW enumerated all codewords to collect (outC, maxC) */
  STOP_DSTOP         /* the lower bound reached the stop target: dmin >= dstop */
};
 
/* Stop all threads, recording the `reason` unless one is already recorded */
static inline void ctx_stop(distfork_ctx_t * const ctx, const int reason) {
  int none = STOP_NONE;
  atomic_compare_exchange_strong(&ctx->stop_reason, &none, reason);
  atomic_store(&ctx->stop_flag, true);
}
 
/* The target t of the CC rounds: the upper bound dmax, or the stop target dstop if smaller (0: neither).  The rounds
 * w <= t-1 certify dmin = t (with outC, the rounds w = t..t+dW collect codewords).  Unlike dmax, dstop is not an upper
 * bound: once dmin >= dstop, the run ends with the bounds [dmin, dmax] (dmax = 0 if no codeword was found). */
static inline int ctx_target(const distfork_ctx_t * const ctx, const int cur_dmax) {
  const int dstop = ctx->p->dstop;
  return (dstop > 0 && (cur_dmax == 0 || dstop < cur_dmax)) ? dstop : cur_dmax;
}
 
typedef struct worker_arg {
  distfork_ctx_t *ctx;
  int tid;
  int min_swei[MAX_W];
  cols_ws_t *irr_ws; /* irred=1: the workspace of the test of the CC (and ksub RW) codewords, see cw_irreducible() */
} worker_arg_t;
 
/* A CC round at weight w hands out the work units u = beg*R .. (end+1)*R-1 one at a time, where beg..end are the start
 * indices of the round (see cc_round_beg(), cc_round_end()) and R = B^L (see cc_round_levels()).  The unit u = i*R + k
 * grows the clusters from the start index i = u / R, with the branches of the first L levels fixed by the digits
 * k_0, k_1, ... of k = k_0 + k_1*B + ... in radix B, the maximum row weight of H: at level l = 1..L (a cluster of l
 * columns), start_CC_recurs_mt() only adds the column at position k_{l-1} of the row it extends (none if the row is
 * shorter: an empty unit).  The units of a start index explore together the same clusters as the start index itself
 * (L = 0), in pieces small enough to balance the load if there are few start indices per thread (the clusters from the
 * first start column, which only grow to larger columns, are about w/n of the round).  The weight, the number of fixed
 * levels L, and the next unit are packed into the single atomic word `cc_next`, so that a worker always processes a
 * claimed unit with the parameters of its round, even if it joins just as the coordinator publishes the next round. */
#define CC_IDX_BITS 32
#define CC_LEV_SHIFT 48
static inline long long cc_pack(const int w, const int lev, const int idx) {
  return ((long long)lev << CC_LEV_SHIFT) | ((long long)w << CC_IDX_BITS) | (long long)(unsigned int)idx;
}
static inline int cc_unpack_w(const long long v) { return (int)((v >> CC_IDX_BITS) & 0xffff); }
static inline int cc_unpack_lev(const long long v) { return (int)((v >> CC_LEV_SHIFT) & 0xf); }
static inline int cc_unpack_idx(const long long v) { return (int)(v & 0xffffffffLL); }
 
/* First start index of a CC round: index into the expert `start` list, or the first start column (cbeg) */
static inline int cc_round_beg(const params_t * const p) {
  return (p->start_num > 0) ? 0 : ((p->cbeg >= 0) ? p->cbeg : 0);
}
 
/* Last start index of a CC round at weight w: a cluster started at column i only uses columns > i, so it
 * starts at column n-w at the latest (expert `start` list: all listed columns, clusters are unlimited) */
static inline int cc_round_end(const params_t * const p, const int w) {
  if (p->start_num > 0) return p->start_num - 1;
  const int nvar = p->spaH->cols;
  return (p->cend >= 0) ? minint(p->cend, nvar - w) : (nvar - w);
}
 
/* Last work unit of the CC round at weight w with `lev` fixed levels (in range, see cc_round_levels()) */
static inline int cc_unit_last(const distfork_ctx_t * const ctx, const int w, const int lev) {
  return (cc_round_end(ctx->p, w) + 1) * ctx->cc_R[lev] - 1;
}
 
#define CC_SPLIT_UNITS_PER_THREAD 4.0 /* aim at >= 4 work units of the size of the largest one per CC thread */
#define CC_SPLIT_MIN_WORK 1e-3        /* rounds predicted to take < 1 ms on all threads (on the start indices of an
                                       * expert `start` list, if they are fewer): whole start indices */
#define CC_SPLIT_CLAIMS_PER_WORK 1e4  /* at most 10^4 work units per thread-second of predicted work */
#define CC_SPLIT_CLAIMS_UNKNOWN 65536 /* at most this many work units if the work of the round is not known */
 
/* Number of fixed levels L of the work units of the CC round at weight w over the start indices beg..end, with the
 * estimated work `est` (thread-seconds, 0 if not known, see cc_work_estimate()), see `cc_pack()`: the smallest L which
 * gives CC_SPLIT_UNITS_PER_THREAD units per thread, where the largest unit is taken as the unit size: whole start
 * indices count as n/w units (the clusters from the first start column are about w/n of the round; the columns of the
 * expert `start` list are alike, with unlimited clusters), and each fixed level divides the units by the typical
 * number of branches.  L = 0 for short rounds (with a short `start` list, only if they are short on its start indices,
 * which are then all that can work in parallel), and L <= min(3, w-3), so that at least two levels of a unit are free.
 * The number of units (each one claim of the shared counter, empty units included) is limited by the predicted work,
 * and all units must have an index < INT_MAX - num_threads (over-claims at the end of the round). */
static int cc_round_levels(const distfork_ctx_t * const ctx, const int w, const int beg, const int end,
                           const double est) {
  const int nthr = ctx->num_threads;
  const int n_start = end - beg + 1;
  const int lmax = minint(CC_SPLIT_MAX_LEVELS, w - 3);
  if (nthr <= 1 || n_start <= 0 || lmax <= 0 || ctx->cc_B <= 1) return 0;
  const int busy = (ctx->p->start_num > 0) ? minint(nthr, n_start) : nthr; /* threads busy with whole start indices */
  if (est > 0.0 && est < CC_SPLIT_MIN_WORK * busy) return 0;
  const double max_units = (est > 0.0) ? CC_SPLIT_CLAIMS_PER_WORK * est : CC_SPLIT_CLAIMS_UNKNOWN;
  double units = (ctx->p->start_num > 0) ? (double)n_start : (double)n_start / w;
  int lev = 0;
  while (lev < lmax && units < CC_SPLIT_UNITS_PER_THREAD * nthr) {
    const long long R = ctx->cc_R[lev + 1];
    if (R <= 0 || (double)n_start * (double)R > max_units || ((long long)end + 1) * R > (long long)INT_MAX - nthr - 1)
      break;
    lev++;
    units *= ctx->cc_branch;
  }
  return lev;
}
 
/* Update the RW weight limit for the `min_hits` statistic after a change of the codeword hash (under `cw_mutex`).  It
 * is stored only if changed, as it shares the cache line of the bounds, which RW threads read for every candidate. */
static inline void ctx_update_feed_limit(distfork_ctx_t * const ctx) {
  const int feed = codeword_feed_limit(ctx->p);
  if (atomic_load_explicit(&ctx->cw_feed_limit, memory_order_relaxed) != feed) {
    atomic_store_explicit(&ctx->cw_feed_limit, feed, memory_order_relaxed);
  }
}
 
/* irred=1: is the codeword of the `cnt` columns `vec` irreducible, i.e., not the union of a codeword and a disjoint
 * nonzero zero-syndrome vector?  Iff the rank of its columns of H is cnt - 1 (see csr_cols_rank()).  Always true with
 * irred=0. */
static inline bool cw_irreducible(const distfork_ctx_t * const ctx, cols_ws_t * const ws, const int cnt,
                                  const int vec[]) {
  if (!ctx->p->irred || !ws || cnt < 2) return true;
  return csr_cols_rank(ctx->mHT_cc, cnt, vec, ws) == cnt - 1;
}
 
/* Record the CC codeword `err` (sorted columns, nonzero syndrome of L) found by worker `warg`, and update the upper
 * bound; returns true if the run is to stop (the first codeword without outC and maxC, or maxC codewords).  With
 * irred=1, a reducible codeword is not recorded (counted in `cw_red_cc`); the CC search finds every irreducible
 * codeword, but only some reducible ones, from the weight d + 3 on for a circuit DEM (doc/cc_completeness.tex). */
static __attribute__((noinline)) bool cc_record(distfork_ctx_t * const ctx, const worker_arg_t * const warg,
                                                const one_vec_t * const err) {
  params_t * const p = ctx->p;
  if (!cw_irreducible(ctx, warg->irr_ws, err->wei, err->vec)) {
    atomic_fetch_add_explicit(&ctx->cw_red_cc, 1, memory_order_relaxed);
    return false;
  }
  bool stop = false;
  pthread_mutex_lock(&ctx->cw_mutex);
  p->codewords = codeword_add_maybe(p, err->vec, err->wei, CW_SRC_OTHER);
  ctx_update_feed_limit(ctx);
  int cur_d = atomic_load(&ctx->dmax);
  if (p->min_w < cur_d || cur_d == 0) {
    atomic_store(&ctx->dmax, p->min_w);
    if ((p->debug & DBG_CODEWORDS) && err->wei == p->min_w) {
      char prefix[48];
      snprintf(prefix, sizeof(prefix), "# [thread %d] CC codeword", warg->tid);
      print_codeword_support(stderr, prefix, err->vec, err->wei, p->debug);
    }
  }
  int cur_cc_found = atomic_load(&ctx->cc_found_weight);
  if (cur_cc_found == 0 || err->wei < cur_cc_found) {
    atomic_store(&ctx->cc_found_weight, err->wei);
  }
  if (!p->outC && p->maxC == 0) {
    ctx_stop(ctx, STOP_CC_FOUND);
    stop = true;
  }
  if (codeword_maxc_reached(p)) {
    ctx_stop(ctx, STOP_MAXC);
    stop = true;
  }
  pthread_mutex_unlock(&ctx->cw_mutex);
  return stop;
}
 
static inline int bit_get(const uint64_t * const b, const int i) { return (int)((b[i >> 6] >> (i & 63)) & 1); }
static inline void bit_flip(uint64_t * const b, const int i) { b[i >> 6] ^= (uint64_t)1 << (i & 63); }
 
/* The first set bit of the bitset b in from..limit-1 (limit if none) */
static inline int bits_first(const uint64_t * const b, const int from, const int limit) {
  if (from >= limit) return limit;
  int k = from >> 6;
  const int k_last = (limit - 1) >> 6;
  uint64_t x = b[k] & (~(uint64_t)0 << (from & 63));
  while (x == 0) {
    if (++k > k_last) return limit;
    x = b[k];
  }
  const int r = (k << 6) + __builtin_ctzll(x);
  return (r < limit) ? r : limit;
}
 
/* The state of the CC walk of a worker (smax = 0, see cc_walk() and, for cc_prune >= 1, cc_bnb()) */
typedef struct {
  distfork_ctx_t *ctx;
  worker_arg_t *warg;
  const csr_t *mH;     /* H (sorted rows) */
  const csr_t *mHT;    /* its transpose: the columns of H (sorted rows) */
  const cc_hash_t *hs; /* the column-syndrome hash */
  uint64_t *S;         /* the syndrome of the cluster: bitset of the rows of H */
  uint64_t *E;         /* the columns of the cluster: bitset */
  int *cols;           /* the columns of the cluster in the order added (cols[0]: the start column) */
  one_vec_t *cw;       /* a codeword found: its sorted columns */
  int w;               /* the number of columns of the cluster */
  int col_min;         /* columns are added only if larger (-1: unlimited clusters with the expert `start` list) */
  int w_limit;         /* the weight of the CC round */
  int max_col_wt;      /* the maximum column weight of H */
  int m;               /* the number of rows of H */
  const int *fix;      /* the fixed levels 1..nfix of the work unit, see `cc_pack()` */
  int nfix;
  long long nodes;     /* calls of the recursion (cc_walk() or cc_bnb()), added to ctx->cc_nodes after each round */
  int count_w;         /* the nodes of weight w >= count_w are counted: a node at a fixed level is visited by all the
                        * work units with the same earlier digits, and counted in the one whose later digits are 0 */
  /* the branch and bound (cc_prune >= 1, see cc_bnb()): a column u is available if u > col_min and u is neither in
   * the cluster S nor in the forbidden set F */
  int prune;           /* cc_prune: 1, or 2 with the coverage and fractional cover bounds */
  int max_row_wt;      /* the maximum row weight of H (at most this many branches per node) */
  uint64_t *X;         /* the columns of S and F: bitset */
  int *acnt;           /* the number of available columns of each row of H */
  int acnt_col;        /* the value of col_min for which acnt is set up (INT_MIN: not set up) */
  int *kfirst;         /* the position in mH of the first column u > col_min of each row (the columns of a row are
                        * sorted, so that the loops over the available columns of a row start there) */
  int *fstk;           /* the forbidden columns in the order forbidden (unwound when the node returns) */
  int fsp;             /* their number */
  int *brl;            /* the branches of the nodes: max_row_wt entries per depth */
  int *ulist;          /* the unsatisfied rows of the current node */
  int *uacnt;          /* their numbers of available columns (acnt), for the packing bound */
  int *ord;            /* the unsatisfied rows ordered by their available columns, for the packing bound */
  int *bcnt;           /* counts of the rows by their available columns (max_row_wt + 2 entries) */
  uint32_t *mark;      /* columns marked in generation `gen` (packing and coverage bounds) */
  uint32_t gen;
  int *cv;             /* cc_prune=2: the unsatisfied rows of each available column (generation `gen`) */
  int *touched;        /* cc_prune=2: the available columns of the unsatisfied rows */
  int *hist;           /* cc_prune=2: histogram of cv (max_col_wt + 1 entries) */
} cc_walk_t;
 
/* Record the codeword: the cluster with the columns c1 (unless negative) and c2, if its syndrome of L is nonzero;
 * returns true if the run is to stop */
static bool cc_walk_found(cc_walk_t * const t, const int c1, const int c2) {
  one_vec_t * const cw = t->cw;
  int n = 0;
  for (int i = 0; i < t->w; i++) cw->vec[n++] = t->cols[i];
  if (c1 >= 0) cw->vec[n++] = c1;
  cw->vec[n++] = c2;
  for (int i = 1; i < n; i++) { /* insertion sort */
    const int x = cw->vec[i];
    int j = i - 1;
    while (j >= 0 && cw->vec[j] > x) {
      cw->vec[j + 1] = cw->vec[j];
      j--;
    }
    cw->vec[j + 1] = x;
  }
  cw->wei = n;
  const colmask_t * const mL = t->ctx->p->maskL; /* NULL for a classical code */
  if (mL && !colmask_syndrome_non_zero(mL, n, cw->vec)) return false;
  return cc_record(t->ctx, t->warg, cw);
}
 
/* Is the syndrome S (bitset of weight sw) the sum h_{c1} + h_{c2} of the columns c1 and c2 of H (rows of mHT)? */
static bool cc_sum2_equal(const uint64_t * const S, const int sw, const csr_t * const mHT, const int c1, const int c2) {
  int a = mHT->p[c1], b = mHT->p[c2], k = 0;
  const int a_end = mHT->p[c1 + 1], b_end = mHT->p[c2 + 1];
  while (a < a_end || b < b_end) {
    int r;
    if (b >= b_end) r = mHT->i[a++];
    else if (a >= a_end) r = mHT->i[b++];
    else if (mHT->i[a] < mHT->i[b]) r = mHT->i[a++];
    else if (mHT->i[b] < mHT->i[a]) r = mHT->i[b++];
    else { /* a common row cancels */
      a++;
      b++;
      continue;
    }
    if (!bit_get(S, r)) return false;
    k++;
  }
  return k == sw;
}
 
/* The last two levels of cc_walk() at the level w = t->w with w + 2 = current_limit (the leaf level w + 1 is not a
 * fixed level of the work unit; `i1_beg`..`i1_end`: the positions in mH of the columns of the row r0 = min(S) to add).
 * There, the recursion adds each column c1 > col_min of r0 not in the cluster with the syndrome s = S + h_{c1} nonzero,
 * and the leaf level adds each column c2 > col_min of the row min(s) not in the cluster with h_{c2} = s, which gives a
 * codeword.  These columns c2 are found by a lookup of the hash h ^ hc[c1] of s instead (cc_hash_t, verified exactly).
 * The pairs (c1, c2) are the same and come in the same order (c2 in ascending order, as in the row min(s)).  As at the
 * entry of the leaf level, a codeword is not recorded once the stop flag is set or the weight limit has dropped (a
 * smaller upper bound).  Returns 1 if the run is to stop. */
static __attribute__((noinline)) int cc_walk_last2(cc_walk_t * const t, const int sw, const uint64_t h,
                                                   const int i1_beg, const int i1_end) {
  distfork_ctx_t * const ctx = t->ctx;
  const cc_hash_t * const hs = t->hs;
  const uint64_t * const hc = hs->hc;
  const uint64_t * const bits = hs->bits;
  const uint64_t bmask = hs->bmask;
  const csr_t * const mH = t->mH;
  const int col_min = t->col_min;
  for (int i1 = i1_beg; i1 < i1_end; i1++) {
    const int c1 = mH->i[i1];
    if (c1 <= col_min) continue;
    const uint64_t key = h ^ hc[c1]; /* the hash of s = S + h_{c1} */
    const uint64_t fb = key & bmask;
    if (!((bits[fb >> 6] >> (fb & 63)) & 1)) continue;
    int j = (int)(key >> hs->tshift);
    while (hs->tab[j].cnt > 0 && hs->tab[j].key != key) j = (j + 1) & hs->tmask;
    if (hs->tab[j].cnt == 0) continue;
    /* columns c2 with the hash of s (s is nonzero: no zero columns in the table) */
    if (bit_get(t->E, c1)) continue;
    const int k_end = hs->tab[j].beg + hs->tab[j].cnt;
    for (int k = hs->tab[j].beg; k < k_end; k++) {
      const int c2 = hs->cols[k];
      if (c2 <= col_min || c2 == c1 || bit_get(t->E, c2)) continue;
      if (!cc_sum2_equal(t->S, sw, t->mHT, c1, c2)) continue; /* a hash collision */
      /* the checks at the entry of the leaf level */
      if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed)) return 0;
      const int cur_dmax = atomic_load_explicit(&ctx->dmax, memory_order_relaxed);
      if (cur_dmax > 0 && ctx->p->dW >= 0 && minint(t->w_limit, cur_dmax + ctx->p->dW) < t->w + 2) return 0;
      if (cc_walk_found(t, c1, c2)) return 1;
    }
  }
  return 0;
}
 
/* The CC recursion of start_CC_recurs_mt() for smax = 0, with the syndrome of the cluster as the bitset t->S (weight
 * sw, smallest row r0, hash h, see cc_hash_t), and its columns as the bitset t->E: a column is added by flipping its
 * rows in S (and back after the recursion) instead of merging sorted lists, and the weight of the new syndrome is
 * known before.  The clusters are the same, in the same order, and so are the codewords found (and their CC hit
 * counts).  The last two levels use hash lookups (cc_walk_last2()).  Returns 1 if the run is to stop. */
static int cc_walk(cc_walk_t * const t, const int sw, const int r0, const uint64_t h) {
  distfork_ctx_t * const ctx = t->ctx;
  if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed)) return 0;
  const params_t * const p = ctx->p;
  const int w = t->w;
  int current_limit = t->w_limit;
  const int cur_dmax = atomic_load_explicit(&ctx->dmax, memory_order_relaxed);
  if (cur_dmax > 0 && p->dW >= 0) current_limit = minint(t->w_limit, cur_dmax + p->dW);
  if (w >= current_limit) return 0;
  if (w >= t->count_w) t->nodes++;
  const csr_t * const mH = t->mH;
  const csr_t * const mHT = t->mHT;
  /* the positions in mH of the columns of the row r0 to add: all of them, or one at a fixed level of the work unit */
  int i1_beg = mH->p[r0], i1_end = mH->p[r0 + 1];
  if (w <= t->nfix) {
    i1_beg += t->fix[w - 1];
    if (i1_beg >= i1_end) return 0;
    i1_end = i1_beg + 1;
  }
  const int col_min = t->col_min;
  uint64_t * const S = t->S;
  uint64_t * const E = t->E;
 
  if (w + 1 == current_limit) { /* leaf level: the columns c with h_c = S */
    for (int i1 = i1_beg; i1 < i1_end; i1++) {
      const int c = mH->i[i1];
      if (c <= col_min) continue;
      const int pb = mHT->p[c];
      if (mHT->p[c + 1] - pb != sw) continue;
      int k = 0;
      while (k < sw && bit_get(S, mHT->i[pb + k])) k++;
      if (k < sw || bit_get(E, c)) continue;
      if (cc_walk_found(t, -1, c)) return 1;
    }
    return 0;
  }
  if (w + 2 == current_limit && w + 1 > t->nfix) return cc_walk_last2(t, sw, h, i1_beg, i1_end);
 
  /* internal level */
  const int max_reach = (current_limit - (w + 1)) * t->max_col_wt;
  const uint64_t * const hc = t->hs->hc;
  for (int i1 = i1_beg; i1 < i1_end; i1++) {
    const int c1 = mH->i[i1];
    if (c1 <= col_min) continue;
    const int pb = mHT->p[c1], pe = mHT->p[c1 + 1];
    if (sw - (pe - pb) > max_reach || bit_get(E, c1)) continue;
    int nset = 0; /* the rows of c1 in S */
    for (int k = pb; k < pe; k++) nset += bit_get(S, mHT->i[k]);
    const int sw1 = sw + (pe - pb) - 2 * nset;
    if (sw1 <= 0 || sw1 > max_reach) continue;
    /* add c1: flip its rows; the smallest row of the new syndrome (r0 is in c1, all the other rows of S are larger) */
    int rmin = t->m;
    for (int k = pb; k < pe; k++) {
      const int r = mHT->i[k];
      if (!bit_get(S, r) && r < rmin) rmin = r;
      bit_flip(S, r);
    }
    const int r1 = (rmin < r0) ? rmin : bits_first(S, r0 + 1, rmin);
    bit_flip(E, c1);
    t->cols[t->w++] = c1;
    const int result = cc_walk(t, sw1, r1, h ^ hc[c1]);
    t->w--;
    bit_flip(E, c1);
    for (int k = pb; k < pe; k++) bit_flip(S, mHT->i[k]);
    if (result == 1) return 1;
  }
  return 0;
}
 
/* The CC work unit from the start column `col` (w_limit > 1, see cc_walk()) */
static void cc_walk_unit(cc_walk_t * const t, const int col, const int w_limit, const int * const fix, const int nfix) {
  const csr_t * const mHT = t->mHT;
  const int pb = mHT->p[col], pe = mHT->p[col + 1];
  for (int k = pb; k < pe; k++) bit_flip(t->S, mHT->i[k]);
  bit_flip(t->E, col);
  t->cols[0] = col;
  t->w = 1;
  t->col_min = (t->ctx->p->start_num > 0 && !t->ctx->p->start_scan) ? -1 : col;
  t->w_limit = w_limit;
  t->fix = fix;
  t->nfix = nfix;
  t->count_w = nfix + 1; /* (see cc_walk_t: count_w - 1 = the first fixed level with only zero digits from there on) */
  while (t->count_w > 1 && fix[t->count_w - 2] == 0) t->count_w--;
  cc_walk(t, pe - pb, mHT->i[pb], t->hs->hc[col]);
  t->w = 0;
  bit_flip(t->E, col);
  for (int k = pb; k < pe; k++) bit_flip(t->S, mHT->i[k]);
}
 
/* ------------------------------------------------------------------------------------------------------------------
 * The branch-and-bound CC search (cc_prune >= 1, smax = 0; see doc/cc_completeness.tex for the proof).  A node is a
 * cluster S with a nonzero syndrome s and a forbidden set F: a column u is available if u > col_min (the start column,
 * F_0 = {u < start column}; with the expert `start` list, the columns listed before the start column) and u is not
 * in S or F.  With R = limit - |S| columns left, a node is pruned by admissible lower bounds on the number of columns
 * needed to complete S to a codeword: |s| > R * max_col_wt (b1, as in cc_walk()), an unsatisfied row without
 * available columns (dead end), more than R unsatisfied rows with pairwise disjoint sets of available columns (b2,
 * packing, found greedily; with R >= 3), and with cc_prune=2 the coverage and fractional cover bounds (with R >= 4).
 * Otherwise it branches on the unsatisfied row with the fewest available columns (the first such row in the row
 * order), over its available columns v_1..v_q in column order: the branch i adds v_i to S, and v_i is forbidden in the
 * later branches whatever happened in its own (sibling exclusion).  Every irreducible codeword which contains the
 * start column and avoids F_0 is found exactly once (one CC hit per codeword), reducible ones at most once.  The
 * choice of the row and the order of the branches depend only on the node, so that the work units (the fixed levels,
 * see `cc_pack()`: the branch fix[l] among the available columns of the row, with the earlier ones forbidden)
 * partition the search.  The last two levels use the column-syndrome hash, as in cc_walk_last2().
 * ------------------------------------------------------------------------------------------------------------------ */
 
/* The column v becomes unavailable (added to S or to F) */
static inline void bnb_take(cc_walk_t * const t, const int v) {
  const csr_t * const mHT = t->mHT;
  const int * const ri = mHT->i, ke = mHT->p[v + 1];
  int * const acnt = t->acnt;
  t->X[v >> 6] |= (uint64_t)1 << (v & 63);
  for (int k = mHT->p[v]; k < ke; k++) acnt[ri[k]]--;
}
 
/* The column v becomes available again */
static inline void bnb_give(cc_walk_t * const t, const int v) {
  const csr_t * const mHT = t->mHT;
  const int * const ri = mHT->i, ke = mHT->p[v + 1];
  int * const acnt = t->acnt;
  t->X[v >> 6] &= ~((uint64_t)1 << (v & 63));
  for (int k = mHT->p[v]; k < ke; k++) acnt[ri[k]]++;
}
 
static inline bool bnb_avail(const cc_walk_t * const t, const int u) { return u > t->col_min && !bit_get(t->X, u); }
 
/* bnb_avail() for a column u > col_min (from the position t->kfirst[r] on in its row r) */
static inline bool bnb_free(const cc_walk_t * const t, const int u) { return !bit_get(t->X, u); }
 
/* A new generation of the column marks (cleared when the counter wraps around) */
static inline uint32_t bnb_next_gen(cc_walk_t * const t) {
  if (++t->gen == 0) {
    memset(t->mark, 0, (size_t)t->mHT->rows * sizeof(uint32_t));
    t->gen = 1;
  }
  return t->gen;
}
 
/* As at the entry of a level: no codeword of weight |S| + extra is recorded once the stop flag is set or the weight
 * limit has dropped below it (a smaller upper bound) */
static inline bool bnb_record_ok(const cc_walk_t * const t, const int extra) {
  distfork_ctx_t * const ctx = t->ctx;
  if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed)) return false;
  const int cur_dmax = atomic_load_explicit(&ctx->dmax, memory_order_relaxed);
  return !(cur_dmax > 0 && ctx->p->dW >= 0 && minint(t->w_limit, cur_dmax + ctx->p->dW) < t->w + extra);
}
 
/* The packing bound b2 for the nu unsatisfied rows t->ulist (t->uacnt: their amin..amax available columns): a greedy
 * packing of rows with pairwise disjoint sets of available columns, the rows with the fewest available columns first
 * (in row order); returns its size once it exceeds R, or once the rows left cannot make it exceed R */
static int bnb_packing(cc_walk_t * const t, const int R, const int nu, const int amin, const int amax) {
  int * const bcnt = t->bcnt; /* bcnt[a - amin]: first the number of the rows with a - 1 available columns */
  const int * const ulist = t->ulist, * const uacnt = t->uacnt;
  int * const ord = t->ord;
  const int na = amax - amin + 1;
  for (int a = 0; a <= na; a++) bcnt[a] = 0;
  for (int i = 0; i < nu; i++) bcnt[uacnt[i] - amin + 1]++;
  for (int a = 1; a <= na; a++) bcnt[a] += bcnt[a - 1]; /* then the first position of the rows with a */
  for (int i = 0; i < nu; i++) ord[bcnt[uacnt[i] - amin]++] = ulist[i];
  const uint32_t g = bnb_next_gen(t);
  uint32_t * const mark = t->mark; /* the columns of the rows packed so far (the available ones count) */
  const uint64_t * const X = t->X;
  const int * const rp = t->mH->p, * const ci = t->mH->i, * const kf = t->kfirst;
  int cnt = 0;
  for (int i = 0; i < nu; i++) {
    const int r = ord[i], kb = kf[r], ke = rp[r + 1];
    int k = kb;
    while (k < ke && !(mark[ci[k]] == g && !bit_get(X, ci[k]))) k++;
    if (k == ke) { /* disjoint from the rows packed so far */
      if (++cnt > R || cnt + (nu - 1 - i) <= R) return cnt;
      for (k = kb; k < ke; k++) mark[ci[k]] = g;
    } else if (cnt + (nu - 1 - i) <= R) {
      return cnt;
    }
  }
  return cnt;
}
 
/* cc_prune=2: the coverage bound b3 (the least t such that the t largest values c_u of the available columns sum to at
 * least |s|, where c_u is the number of unsatisfied rows of u) and the fractional cover bound b4 (the sum over the
 * unsatisfied rows r of 1 / max c_u over the available columns u of r); true if one of them exceeds R */
static bool bnb_cover(cc_walk_t * const t, const int R, const int nu) {
  const csr_t * const mH = t->mH;
  const uint32_t g = bnb_next_gen(t);
  int nt = 0;
  for (int i = 0; i < nu; i++) {
    const int r = t->ulist[i];
    for (int k = t->kfirst[r]; k < mH->p[r + 1]; k++) {
      const int u = mH->i[k];
      if (!bnb_free(t, u)) continue;
      if (t->mark[u] != g) {
        t->mark[u] = g;
        t->cv[u] = 0;
        t->touched[nt++] = u;
      }
      t->cv[u]++;
    }
  }
  const int D = t->max_col_wt;
  int * const hist = t->hist;
  for (int c = 0; c <= D; c++) hist[c] = 0;
  for (int i = 0; i < nt; i++) hist[t->cv[t->touched[i]]]++;
  int need = nu, tc = 0;
  for (int c = D; c >= 1 && need > 0; c--) {
    if ((long long)hist[c] * c >= need) {
      tc += (need + c - 1) / c;
      need = 0;
    } else {
      tc += hist[c];
      need -= hist[c] * c;
    }
  }
  if (need > 0 || tc > R) return true;
  double sum = 0.0; /* (the rows in row order: the same sum on every thread) */
  for (int i = 0; i < nu; i++) {
    const int r = t->ulist[i];
    int mx = 0;
    for (int k = t->kfirst[r]; k < mH->p[r + 1]; k++) {
      const int u = mH->i[k];
      if (bnb_free(t, u) && t->cv[u] > mx) mx = t->cv[u];
    }
    if (mx == 0) return true;
    sum += 1.0 / mx;
    /* the partial sums grow, by at most 1 per row, so that the decision may be known early (the rounding errors are
     * far smaller than the margins) */
    if (sum > R + 1e-9) return true;
    if (sum + (nu - 1 - i) < R + 1e-9 - 1e-7) return false;
  }
  return false;
}
 
/* The hash slot of the key (a syndrome hash, see cc_hash_t), NULL if no nonzero column has this hash */
static inline const cc_hash_slot_t *cc_hash_find(const cc_hash_t * const hs, const uint64_t key) {
  const uint64_t fb = key & hs->bmask;
  if (!((hs->bits[fb >> 6] >> (fb & 63)) & 1)) return NULL;
  int j = (int)(key >> hs->tshift);
  while (hs->tab[j].cnt > 0 && hs->tab[j].key != key) j = (j + 1) & hs->tmask;
  return (hs->tab[j].cnt > 0) ? &hs->tab[j] : NULL;
}
 
/* The leaf level of cc_bnb() (one column left, not a fixed level): the available columns u with h_u = s, found by a
 * lookup of the hash h of s (verified exactly) */
static int bnb_leaf(cc_walk_t * const t, const int sw, const uint64_t h) {
  const cc_hash_slot_t * const sl = cc_hash_find(t->hs, h);
  if (!sl) return 0;
  const csr_t * const mHT = t->mHT;
  for (int k = sl->beg; k < sl->beg + sl->cnt; k++) {
    const int u = t->hs->cols[k];
    if (!bnb_avail(t, u)) continue;
    const int pb = mHT->p[u];
    if (mHT->p[u + 1] - pb != sw) continue;
    int x = 0;
    while (x < sw && bit_get(t->S, mHT->i[pb + x])) x++;
    if (x < sw) continue;
    if (!bnb_record_ok(t, 1)) return 0;
    if (cc_walk_found(t, -1, u)) return 1;
  }
  return 0;
}
 
/* The last two levels of cc_bnb() (two columns left, not a fixed level): for each available column v of the branching
 * row r (in column order), the available columns u != v with h_u = s + h_v, found by a lookup of the hash h ^ hc[v]
 * (verified exactly), give the codewords S + {v, u}, each exactly once (doc/cc_completeness.tex, Proposition "the last
 * two extensions by a hash table"; the columns v of the earlier branches need not be forbidden: such a column u has
 * a zero in the row r, unlike v) */
static int bnb_last2(cc_walk_t * const t, const int sw, const uint64_t h, const int r) {
  const cc_hash_t * const hs = t->hs;
  const csr_t * const mH = t->mH;
  for (int i = t->kfirst[r]; i < mH->p[r + 1]; i++) {
    const int v = mH->i[i];
    if (!bnb_free(t, v)) continue;
    const cc_hash_slot_t * const sl = cc_hash_find(hs, h ^ hs->hc[v]);
    if (!sl) continue;
    for (int k = sl->beg; k < sl->beg + sl->cnt; k++) {
      const int u = hs->cols[k];
      if (u == v || !bnb_avail(t, u)) continue;
      if (!cc_sum2_equal(t->S, sw, t->mHT, v, u)) continue; /* a hash collision */
      if (!bnb_record_ok(t, 2)) return 0;
      if (cc_walk_found(t, v, u)) return 1;
    }
  }
  return 0;
}
 
/* The branch-and-bound recursion at the node (S, F) with S = t->cols[0..w-1] (syndrome t->S of weight sw > 0, hash
 * h), F = the columns of t->X not in S; returns 1 if the run is to stop.  On return, t->X, t->acnt, and t->fsp are as
 * on entry. */
static int cc_bnb(cc_walk_t * const t, const int sw, const uint64_t h) {
  distfork_ctx_t * const ctx = t->ctx;
  if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed)) return 0;
  const params_t * const p = ctx->p;
  const int w = t->w;
  int current_limit = t->w_limit;
  const int cur_dmax = atomic_load_explicit(&ctx->dmax, memory_order_relaxed);
  if (cur_dmax > 0 && p->dW >= 0) current_limit = minint(t->w_limit, cur_dmax + p->dW);
  const int R = current_limit - w; /* the number of columns left */
  if (R <= 0 || sw > R * t->max_col_wt) return 0; /* b1 */
  if (w >= t->count_w) t->nodes++;
  const bool fixed = (w <= t->nfix); /* a fixed level of the work unit */
  if (R == 1 && !fixed) return bnb_leaf(t, sw, h);
  /* the unsatisfied rows, and the branching row: the one with the fewest available columns, the first in row order */
  const int nwords = (t->m + 63) >> 6;
  const int * const acnt = t->acnt;
  int nu = 0, rbest = -1, abest = INT_MAX, amax = 0;
  if (R == 2 && !fixed) { /* (only the branching row is needed) */
    for (int k = 0; k < nwords; k++) {
      uint64_t x = t->S[k];
      while (x) {
        const int r = (k << 6) + __builtin_ctzll(x);
        x &= x - 1;
        const int a = acnt[r];
        if (a < abest) {
          if (a == 0) return 0; /* a dead end */
          abest = a;
          rbest = r;
        }
      }
    }
    return bnb_last2(t, sw, h, rbest);
  }
  int * const ulist = t->ulist, * const uacnt = t->uacnt;
  for (int k = 0; k < nwords; k++) {
    uint64_t x = t->S[k];
    while (x) {
      const int r = (k << 6) + __builtin_ctzll(x);
      x &= x - 1;
      const int a = acnt[r];
      ulist[nu] = r;
      uacnt[nu++] = a;
      if (a < abest) {
        abest = a;
        rbest = r;
      }
      if (a > amax) amax = a;
    }
  }
  if (abest == 0) return 0; /* a dead end */
  if (sw > R && R >= 3) { /* (the packing has at most |s| rows; with R <= 2, the hash levels are cheaper) */
    if (bnb_packing(t, R, nu, abest, amax) > R) return 0;
    /* (with R = 3, the cover bounds rarely prune what the packing leaves, and cost more than they save) */
    if (t->prune >= 2 && R >= 4 && bnb_cover(t, R, nu)) return 0;
  }
  /* the branches: the available columns of the row rbest, in column order */
  const csr_t * const mH = t->mH;
  const csr_t * const mHT = t->mHT;
  int * const brl = t->brl + (size_t)w * (size_t)t->max_row_wt;
  int q = 0;
  for (int k = t->kfirst[rbest]; k < mH->p[rbest + 1]; k++) {
    const int u = mH->i[k];
    if (bnb_free(t, u)) brl[q++] = u;
  }
  const int fsp0 = t->fsp;
  int i_beg = 0, i_end = q;
  if (fixed) { /* the branch fix[w-1] only, with the columns of the earlier branches forbidden */
    const int d = t->fix[w - 1];
    if (d >= q) return 0;
    for (int i = 0; i < d; i++) {
      bnb_take(t, brl[i]);
      t->fstk[t->fsp++] = brl[i];
    }
    i_beg = d;
    i_end = d + 1;
  }
  const int reach = (R - 1) * t->max_col_wt; /* the largest syndrome weight which R - 1 columns can complete */
  const uint64_t * const hc = t->hs->hc;
  uint64_t * const S = t->S;
  int res = 0;
  for (int i = i_beg; i < i_end && !res; i++) {
    const int v = brl[i];
    const int pb = mHT->p[v], pe = mHT->p[v + 1];
    bnb_take(t, v); /* v is in S in its branch, and in F after it, whatever happens in the branch */
    t->fstk[t->fsp++] = v;
    if (sw - (pe - pb) > reach) continue; /* b1 of the branch: its syndrome weight is at least sw - |h_v| */
    int nset = 0; /* the rows of v in s */
    for (int k = pb; k < pe; k++) nset += bit_get(S, mHT->i[k]);
    const int sw1 = sw + (pe - pb) - 2 * nset;
    if (sw1 == 0) { /* a codeword if the cluster has the weight of the limit; nothing is searched below it */
      if (R == 1 && bnb_record_ok(t, 1)) res = cc_walk_found(t, -1, v);
      continue;
    }
    if (sw1 > reach) continue;
    for (int k = pb; k < pe; k++) bit_flip(S, mHT->i[k]);
    t->cols[t->w++] = v;
    res = cc_bnb(t, sw1, h ^ hc[v]);
    t->w--;
    for (int k = pb; k < pe; k++) bit_flip(S, mHT->i[k]);
  }
  while (t->fsp > fsp0) bnb_give(t, t->fstk[--t->fsp]);
  return res;
}
 
/* Set up t->acnt and t->kfirst for col_min = c (S and F empty): the columns u > c of each row.  The start columns
 * claimed by a worker grow within a round, so that both are updated incrementally from a smaller c (the columns
 * c' < u <= c are dropped: they are the first columns > c' of their rows); otherwise (a new round, or the expert
 * `start` list after other units) they are counted again. */
static void bnb_acnt_setup(cc_walk_t * const t, const int c) {
  if (t->acnt_col == c) return;
  if (t->acnt_col != INT_MIN && t->acnt_col < c) {
    const csr_t * const mHT = t->mHT;
    for (int u = t->acnt_col + 1; u <= c; u++)
      for (int k = mHT->p[u]; k < mHT->p[u + 1]; k++) {
        t->acnt[mHT->i[k]]--;
        t->kfirst[mHT->i[k]]++;
      }
  } else {
    const csr_t * const mH = t->mH;
    for (int r = 0; r < t->m; r++) {
      int k = mH->p[r];
      while (k < mH->p[r + 1] && mH->i[k] <= c) k++; /* (the columns of a row are sorted) */
      t->acnt[r] = mH->p[r + 1] - k;
      t->kfirst[r] = k;
    }
  }
  t->acnt_col = c;
}
 
/* The CC work unit from the start column `col` with the branch and bound (cc_prune >= 1, w_limit > 1, see cc_bnb()):
 * the clusters use only the columns u > col, or with the expert `start` list (unlimited clusters), all columns except
 * those listed before `col` (`idx`: its position in the list), which are forbidden: a codeword which contains one of
 * them is found from the first one it contains (doc/cc_completeness.tex, Proposition "start lists").  With
 * start_scan=1, the listed start columns use only the larger columns, as in the scan of all columns. */
static void cc_bnb_unit(cc_walk_t * const t, const int col, const int idx, const int w_limit, const int * const fix,
                        const int nfix) {
  const params_t * const p = t->ctx->p;
  const csr_t * const mHT = t->mHT;
  const bool list = (p->start_num > 0 && !p->start_scan);
  t->col_min = list ? -1 : col;
  bnb_acnt_setup(t, t->col_min);
  t->fsp = 0;
  if (list) {
    for (int i = 0; i < idx; i++) {
      bnb_take(t, p->start_list[i]);
      t->fstk[t->fsp++] = p->start_list[i];
    }
    bnb_take(t, col); /* (with col_min = col, the start column is not available anyway) */
  }
  const int pb = mHT->p[col], pe = mHT->p[col + 1];
  for (int k = pb; k < pe; k++) bit_flip(t->S, mHT->i[k]);
  t->cols[0] = col;
  t->w = 1;
  t->w_limit = w_limit;
  t->fix = fix;
  t->nfix = nfix;
  t->count_w = nfix + 1; /* (see cc_walk_t: count_w - 1 = the first fixed level with only zero digits from there on) */
  while (t->count_w > 1 && fix[t->count_w - 2] == 0) t->count_w--;
  cc_bnb(t, pe - pb, t->hs->hc[col]);
  t->w = 0;
  for (int k = pb; k < pe; k++) bit_flip(t->S, mHT->i[k]);
  if (list) bnb_give(t, col);
  while (t->fsp > 0) bnb_give(t, t->fstk[--t->fsp]);
}
 
/* Recursive CC worker function (interruptible), used with smax > 0 (otherwise, cc_walk() explores the same clusters).
 * At the fixed levels w = err->wei <= nfix of a CC work unit (see `cc_pack()`), only the column at position fix[w-1]
 * of the row is added (none if the row is shorter). */
static int start_CC_recurs_mt(one_vec_t *err, one_vec_t *urr, one_vec_t * const syn[],
                              const int w_limit, const int max_col_wt,
                              const csr_t * const mH, const csr_t * const mHT,
                              worker_arg_t *warg, const int * const fix, const int nfix) {
  distfork_ctx_t *ctx = warg->ctx;
  if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed)) {
    return 0;
  }
  params_t * const p = ctx->p;
  const int w = err->wei;
  int current_limit = w_limit;
  int cur_dmax = atomic_load_explicit(&ctx->dmax, memory_order_relaxed);
  if (cur_dmax > 0 && p->dW >= 0) {
    current_limit = minint(w_limit, cur_dmax + p->dW);
  }
  if (w >= current_limit) {
    return 0;
  }
 
  const one_vec_t * const syn_w = syn[w];
  const int syn_w_wei = syn_w->wei;
  const int row = syn_w->vec[0];
  /* the positions in mH of the columns of the row to add: all of them, or one at a fixed level of the work unit */
  int i1_beg = mH->p[row], i1_end = mH->p[row + 1];
  if (w <= nfix) {
    i1_beg += fix[w - 1];
    if (i1_beg >= i1_end) return 0;
    i1_end = i1_beg + 1;
  }
  const colmask_t * const mL = p->maskL; /* NULL for a classical code */
  /* Clusters are grown only to columns larger than the start column, except with the expert `start`
   * list, where clusters are unlimited (columns already in `err` are skipped via `one_ordered_search`),
   * unless start_scan=1. */
  const int col_min = (p->start_num > 0 && !p->start_scan) ? -1 : urr->vec[0];
 
  /* Leaf level: w + 1 == current_limit */
  if (w + 1 == current_limit) {
    int max_leaf_swei = 0;
    if (p->smax > 0 && current_limit < MAX_W) {
      int cur_min = warg->min_swei[current_limit] - 1;
      max_leaf_swei = (cur_min < p->smax) ? (cur_min > 0 ? cur_min : 0) : p->smax;
    }
 
    for (int i1 = i1_beg; i1 < i1_end; i1++) {
      const int col = mH->i[i1];
      if (col <= col_min) continue;
 
      const int p_beg = mHT->p[col];
      const int col_wt = mHT->p[col + 1] - p_beg;
      if (abs(syn_w_wei - col_wt) > max_leaf_swei) continue;
 
      int swei;
      if (max_leaf_swei == 0) {
        /* Only swei == 0 is of interest: check if mHT[col] == syn[w] */
        if (mHT->i[p_beg] != row ||
            mHT->i[p_beg + col_wt - 1] != syn_w->vec[col_wt - 1]) {
          continue;
        }
        if (memcmp(syn_w->vec, &mHT->i[p_beg], (size_t)col_wt * sizeof(int)) != 0) {
          continue;
        }
        if (one_ordered_search(err, col) != -1) continue;
        swei = 0;
      } else {
        if (one_ordered_search(err, col) != -1) continue;
        syn[w + 1]->wei = 0;
        swei = one_csr_row_combine(syn[w + 1], syn_w, mHT, col);
        if (swei > 0) {
          if (swei <= max_leaf_swei) {
            warg->min_swei[w + 1] = swei;
            int cur_min = swei - 1;
            max_leaf_swei = (cur_min < p->smax) ? (cur_min > 0 ? cur_min : 0) : p->smax;
          }
          continue;
        }
      }
 
      /* swei == 0: insert col into err to verify against mL and record codeword */
      int pos = one_ordered_ins(err, col);
      int nz = (!mL) || colmask_syndrome_non_zero(mL, err->wei, err->vec);
      if (nz && cc_record(ctx, warg, err)) {
        one_ordered_pos_del(err, col, pos);
        return 1;
      }
      one_ordered_pos_del(err, col, pos);
    }
    return 0;
  }
 
  /* Internal level: w + 1 < current_limit */
  const int rem = current_limit - (w + 1);
  int max_s_needed = 0;
  if (p->smax > 0 && current_limit < MAX_W) { /* smax > 0 only without noscan and with dmin <= 1 (var_init) */
    int cur_min = warg->min_swei[current_limit] - 1;
    max_s_needed = (cur_min < p->smax) ? (cur_min > 0 ? cur_min : 0) : p->smax;
  }
  const int max_reach = rem * max_col_wt + max_s_needed;
 
  for (int i1 = i1_beg; i1 < i1_end; i1++) {
    const int col = mH->i[i1];
    if (col > col_min) {
      const int col_wt = mHT->p[col + 1] - mHT->p[col];
      if (syn_w_wei - col_wt > max_reach) continue;
 
      int pos = one_ordered_search(err, col);
      if (pos == -1) {
        syn[w + 1]->wei = 0;
        int swei = one_csr_row_combine(syn[w + 1], syn_w, mHT, col);
 
        if (p->smax && swei > 0 && swei <= p->smax && (w + 1 < MAX_W)) {
          if (swei < warg->min_swei[w + 1]) {
            warg->min_swei[w + 1] = swei;
          }
        }
 
        if (swei > 0 && swei <= max_reach) {
          urr->vec[w] = col;
          urr->wei++;
          pos = one_ordered_ins(err, col);
          int result = start_CC_recurs_mt(err, urr, syn, w_limit, max_col_wt,
                                          mH, mHT, warg, fix, nfix);
          urr->wei--;
          one_ordered_pos_del(err, col, pos);
          if (result == 1) {
            return 1;
          }
        }
      }
    }
  }
  return 0;
}
 
/* The window of a hybrid RW step (kwin=0, n >= 500: an index window on every other global step), see
 * rw_hybrid_window() in util_io.h */
static inline int hybrid_window(const params_t * const p, const int nvar) {
  return rw_hybrid_window(p, nvar);
}
 
/* Add the thread time `*acc_t` of `*acc_n` completed RW steps to the totals of the RW timing model, and reset them */
static inline void rw_time_publish(distfork_ctx_t *ctx, double *acc_t, int *acc_n) {
  if (*acc_n > 0) {
    atomic_fetch_add_explicit(&ctx->rw_busy_ns, (long long)(*acc_t * 1e9), memory_order_relaxed);
    atomic_fetch_add_explicit(&ctx->rw_timed_steps, *acc_n, memory_order_relaxed);
  }
  *acc_t = 0.0;
  *acc_n = 0;
}
 
/* Time of a completed RW step since `*t_last`: published at least every 2 ms (after every step if steps are slow),
 * and right away for the first step of a batch (`first`), so that the coordinator gets the step time early */
static inline void rw_time_step(distfork_ctx_t *ctx, double *t_last, double *acc_t, int *acc_n, const bool first) {
  const double t_now = get_time_sec();
  *acc_t += t_now - *t_last;
  *t_last = t_now;
  (*acc_n)++;
  if (first || *acc_t >= 0.002) rw_time_publish(ctx, acc_t, acc_n);
}
 
/* Weight limit (exclusive) of the RW candidate codewords, see codeword_rw_limit() */
static inline int rw_cw_limit(distfork_ctx_t * const ctx) {
  return codeword_rw_limit(ctx->p, atomic_load_explicit(&ctx->dmax, memory_order_relaxed),
                           atomic_load_explicit(&ctx->cw_feed_limit, memory_order_relaxed));
}
 
/* Record the RW codeword `ee` of weight `cnt` (sorted indices) found in a step with a uniform random permutation or
 * with a window (`src`: CW_SRC_RW_UNIFORM or CW_SRC_RW_WINDOW): add it to the hash (or count one more hit), update
 * the upper bound dmax and the stop flags (`ksub` > 0 only for the message) */
static void rw_record_codeword(distfork_ctx_t * const ctx, const rci_t * const ee, const int cnt, const int tid,
                               const int ksub, const int src) {
  params_t * const p = ctx->p;
  /* with outC or maxC, codewords are collected: RW does not end the run when the bounds coincide or on wmin (in
   * method 3, the coordinator then runs the CC rounds which enumerate the codewords) */
  const bool collecting = (p->outC != NULL) || (p->maxC > 0);
  pthread_mutex_lock(&ctx->cw_mutex);
  p->codewords = codeword_add_maybe(p, ee, cnt, src);
  ctx_update_feed_limit(ctx);
  const int best = p->min_w; /* a codeword lighter than min_w is always added */
  const int old_dmax = atomic_load(&ctx->dmax);
  if (old_dmax == 0 || best < old_dmax) {
    atomic_store(&ctx->dmax, best);
    if (p->debug & DBG_PROGRESS) {
      int num_rw = (p->method == 1) ? ctx->num_threads
                   : minint(ctx->rw_threads, ctx->num_threads - atomic_load(&ctx->cc_target_workers));
      if (num_rw < 1) num_rw = 1;
      if (ksub > 0) {
        fprintf(stderr, "# [thread %d] RW (ksub=%d) found new upper bound cw of weight %d (using %d RW threads)\n",
                tid, ksub, best, num_rw);
      } else {
        fprintf(stderr, "# [thread %d] RW found new upper bound cw of weight %d (using %d RW threads)\n",
                tid, best, num_rw);
      }
    }
    if ((p->debug & DBG_CODEWORDS) && cnt == best) {
      char prefix[48];
      snprintf(prefix, sizeof(prefix), "# [thread %d] RW codeword", tid);
      print_codeword_support(stderr, prefix, ee, cnt, p->debug);
    }
    const int cur_dmin = atomic_load(&ctx->dmin);
    if (!collecting && cur_dmin > 0 && best <= cur_dmin) {
      ctx_stop(ctx, STOP_BOUNDS);
    }
  }
  if (!collecting && p->wmin > 0 && best <= p->wmin) {
    ctx_stop(ctx, STOP_WMIN);
  }
  if (codeword_maxc_reached(p)) {
    ctx_stop(ctx, STOP_MAXC);
  }
  if (check_min_hits_convergence(p)) {
    if ((p->debug & DBG_PROGRESS) && !atomic_load(&ctx->stop_flag) && !atomic_load(&ctx->rw_stop_flag)) {
      cw_hit_stats_t st;
      codeword_hit_stats(p, &st);
      fprintf(stderr,
              "# RW convergence reached: <n>=%.2f >= min_hits=%d average hits per codeword (%lld cws of weight "
              "%d..%d; w=%d: %lld cws, %lld hits), P_fail ~ exp(-<n>) = %.2g\n",
              st.avg, p->min_hits, st.set_cws, st.w_lo, st.w_hi, st.w_lo, st.min_cws, st.min_hits, exp(-st.avg));
    }
    if (p->method == 1) {
      ctx_stop(ctx, STOP_RW_CONVERGED);
    } else {
      atomic_store(&ctx->rw_stop_flag, true);
    }
  }
  pthread_mutex_unlock(&ctx->cw_mutex);
}
 
/* dst ^= src, and dst = a ^ b, for rows of `width` words */
static inline void rw_row_xor(word * __restrict dst, const word * __restrict src, const wi_t width) {
  for (wi_t w = 0; w < width; w++) dst[w] ^= src[w];
}
static inline void rw_row_xor3(word * __restrict dst, const word * __restrict a, const word * __restrict b,
                               const wi_t width) {
  for (wi_t w = 0; w < width; w++) dst[w] = a[w] ^ b[w];
}
 
#define RW_K_MAX 8   /* at most 2^8 rows per table of rw_eliminate() */
#define RW_NT_MAX 8  /* at most 8 tables, i.e., 64 pivots per pass over the rows */
 
/* Parameters of rw_eliminate() for a matrix of `rows` x `width` words: the number of tables `*nt` and of pivots per
 * table `*k` (column by column: *nt = 0, below 128 rows).  For wide rows, the tables are kept at about 2 MB or less
 * (cache): fewer rows per table, more tables.  Measured (single thread): K = 6 with 4 tables is 1.7x to 2.2x faster
 * than column by column for 648 to 2736 rows (circuit DEMs, a 1500 x 3000 LDPC code), also with 16 threads; between
 * 128 and 512 rows, K = 5 with 4 tables is 0.75x to 1.4x as fast, depending on the matrix (see rw_elim_t). */
static void rw_block_params(const int rows, const wi_t width, int * const k, int * const nt) {
  *k = 0;
  *nt = 0;
  if (rows < 128) return;
  *k = (rows < 512) ? 5 : 6;
  *nt = 4;
  while (*k > 3 && ((size_t)*nt << *k) * (size_t)width * sizeof(word) > ((size_t)2 << 20)) {
    (*k)--;
    if (*nt + 2 <= RW_NT_MAX) *nt += 2;
  }
}
 
/* RW elimination of a thread: the tables and parameters of rw_eliminate(), and the choice between it and the column-by-
 * column elimination for 128 to 512 rows, where the blocked one is not always faster: the steps 1..2*RW_TUNE_STEPS
 * (after the first one, from the sparse H) alternate between the two eliminations, which give the same results, and
 * the faster one (in total) is used for the following steps */
#define RW_TUNE_STEPS 4
typedef struct {
  word *table;    /* the tables of rw_eliminate(): nt << k rows of tstride words, row 0 of each table zero */
  size_t tstride;
  int k, nt;      /* pivots per table, number of tables (0: column by column), see rw_block_params() */
  int steps;      /* RW steps completed while tuning */
  double t[2];    /* time of the timed eliminations: column by column (0), blocked (1) */
  int choice;     /* -1: tuning, 0: column by column, 1: blocked */
} rw_elim_t;
 
/* RW elimination with the "Method of the Four Russians" (as in M4RI's echelonization, but for the scattered columns
 * of a random permutation): Gauss-Jordan elimination of the columns order[0], order[1], ... in the rows
 * 0..nrows_eff-1 of M, until nrows_eff pivots are found, in blocks of up to nt*K pivots.
 * 1. The pivots of a block, in the order of the columns: a column is a pivot if a row r >= rank, reduced with the
 *    previous pivots, has a 1 in it.  The block's pivot rows are kept reduced with each other, so that the reduced bit
 *    of a row is its bit plus the parity of its bits in the pivot columns of the block's pivot rows with a 1 in the
 *    column: the rows are searched without row operations, and only the new pivot row is reduced.
 * 2. Tables of the 2^K combinations of K pivot rows each (built with one row XOR each, in Gray-code order; row 0 of
 *    each table is zero; table t at T + (t << K) * tstride, see rw_elim_t).
 * 3. Every other row is reduced in one pass: for each table, the combination given by the row's bits in its pivot
 *    columns (the row stays in the L1 cache, while it is read and written only once from memory).
 * The pivots, their rows, and the reduced matrix are the same as with gauss_one_rows() column by column, with about
 * 1/(nt*K) of the passes over the rows of M.  The pivot columns are stored in pivs[0..rank-1] and set in piv_mask.
 * Returns the rank, or -1 if interrupted by a stop flag (checked every 16 columns, as before).  (Not inlined, so that
 * the column-by-column elimination in run_rw_steps() stays as compact as before for small matrices.) */
__attribute__((noinline))
static int rw_eliminate(distfork_ctx_t * const ctx, mzd_t * const M, const int nrows_eff, const int nvar,
                        const rci_t * const order, rci_t * const pivs, word * const piv_mask,
                        const rw_elim_t * const el) {
  word * const T = el->table;
  const size_t tstride = el->tstride;
  const int K = el->k, nt = el->nt;
  const wi_t width = M->width;
  const int P = K * nt; /* pivots per block */
  int bw[RW_K_MAX * RW_NT_MAX], bs[RW_K_MAX * RW_NT_MAX]; /* word and bit of the pivot columns of the block */
  int rank = 0;
  int i = 0;
  while (i < nvar && rank < nrows_eff) {
    const int r0 = rank; /* the first pivot row of the block */
    int kk = 0;
    /* 1. up to P pivots */
    while (kk < P && i < nvar && rank < nrows_eff) {
      if ((i & 15) == 15 && (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed) ||
                             atomic_load_explicit(&ctx->rw_stop_flag, memory_order_relaxed))) {
        return -1;
      }
      const int col = order[i++];
      const int cw = col >> 6, cb = col & 63;
      uint64_t mask = 0; /* the block's pivot rows with a 1 in the column */
      for (int l = 0; l < kk; l++) mask |= (uint64_t)((mzd_row(M, r0 + l)[cw] >> cb) & 1) << l;
      int found = -1;
      for (int r = rank; r < nrows_eff; r++) {
        const word * const row = mzd_row(M, r);
        word b = (row[cw] >> cb) & 1;
        for (uint64_t m2 = mask; m2; m2 &= m2 - 1) {
          const int l = __builtin_ctzll(m2);
          b ^= (row[bw[l]] >> bs[l]) & 1;
        }
        if (b) {
          found = r;
          break;
        }
      }
      if (found < 0) continue; /* not a pivot: the column depends on the previous ones */
      if (found != rank) mzd_row_swap(M, rank, found);
      word * const prow = mzd_row(M, rank);
      for (int l = 0; l < kk; l++) { /* reduce the new pivot row with the block's pivot rows */
        if ((prow[bw[l]] >> bs[l]) & 1) rw_row_xor(prow, mzd_row(M, r0 + l), width);
      }
      for (uint64_t m2 = mask; m2; m2 &= m2 - 1) /* and clear the column in them */
        rw_row_xor(mzd_row(M, r0 + __builtin_ctzll(m2)), prow, width);
      bw[kk] = cw;
      bs[kk] = cb;
      kk++;
      pivs[rank] = col;
      piv_mask[cw] |= (word)1 << cb;
      rank++;
    }
    if (kk == 0) break; /* no more columns */
    /* 2. the tables: row g of table t is the sum of the pivot rows r0 + t*K + l of the bits l of g */
    const int ntab = (kk + K - 1) / K;
    for (int t = 0; t < ntab; t++) {
      word * const Tt = T + ((size_t)t << K) * tstride;
      const int kt = (kk - t * K < K) ? kk - t * K : K;
      for (unsigned g = 1; g < (1u << kt); g++) {
        const unsigned gc = g ^ (g >> 1), gp = (g - 1) ^ ((g - 1) >> 1); /* Gray codes of g and g-1 */
        rw_row_xor3(Tt + gc * tstride, Tt + gp * tstride, mzd_row(M, r0 + t * K + __builtin_ctz(g)), width);
      }
    }
    /* 3. all other rows, in one pass */
    for (int r = 0; r < nrows_eff; r++) {
      if (r == r0) {
        r = rank - 1; /* (the pivot rows of the block) */
        continue;
      }
      word * const row = mzd_row(M, r);
      unsigned x[RW_NT_MAX];
      for (int t = 0; t < ntab; t++) { /* (the bits first: the XORs change the row's bits in the pivot columns) */
        const int lo = t * K, hi = (lo + K < kk) ? lo + K : kk;
        unsigned xt = 0;
        for (int l = lo; l < hi; l++) xt |= (unsigned)((row[bw[l]] >> bs[l]) & 1) << (l - lo);
        x[t] = xt;
      }
      for (int t = 0; t < ntab; t++) {
        if (x[t]) rw_row_xor(row, T + (((size_t)t << K) + x[t]) * tstride, width);
      }
    }
  }
  return rank;
}
 
/* The RW pairs (p=2), see rw_pairs_bits(): a RW step also finds the codewords with exactly two non-pivot columns a, b,
 * the sums of the codewords of a and b (p=1), whose pivot part is the sum of the columns a and b of the reduced matrix.
 * The non-pivot columns are sorted by their bits in l random pivot rows (the bucket, as in Stern's algorithm), and only
 * the pairs in a bucket are checked.  In a bucket, the columns are sorted by weight, so that only the pairs whose
 * weights differ by less than the weight limit are checked (the weight of the sum is at least the difference), and the
 * columns whose codewords have a zero L syndrome come first: the pairs of two of them are skipped (the codewords of
 * most columns of a circuit DEM are light and local, with a zero L syndrome), and so are the pairs with equal L
 * syndromes.  With at most 64 rows of L, the L syndromes of all columns are known (the rows of L reduced with the pivot
 * rows); otherwise, only those of the light columns (a hash, see rw_lhash_init()). */
typedef struct {
  const word *row; /* the column in the pivot rows (a row of the transposed reduced matrix) */
  uint64_t sig;    /* the L syndrome of the column's codeword (or its hash, see rw_lhash_init()), if known */
  int col;         /* the column */
  int known;       /* 1 if `sig` is known */
} rw_pair_t;
 
/* LSD radix sort (8-bit digits) of `n` keys `a` by their lowest `bits` bits (`tmp`: n more) */
static void rw_radix_sort(uint64_t * const a, uint64_t * const tmp, const int n, const int bits) {
  uint64_t *src = a, *dst = tmp;
  for (int sh = 0; sh < bits; sh += 8) {
    int cnt[256];
    memset(cnt, 0, sizeof(cnt));
    for (int i = 0; i < n; i++) cnt[(src[i] >> sh) & 255]++;
    for (int d = 0, sum = 0; d < 256; d++) {
      const int c = cnt[d];
      cnt[d] = sum;
      sum += c;
    }
    for (int i = 0; i < n; i++) dst[cnt[(src[i] >> sh) & 255]++] = src[i];
    uint64_t * const t = src;
    src = dst;
    dst = t;
  }
  if (src != a) memcpy(a, src, (size_t)n * sizeof(uint64_t));
}
 
/* The workspace of a thread for the RW pairs (see rw_pair_t, rw_pairs_step()) */
typedef struct {
  rw_pair_t *ent; /* the entries of the non-pivot columns (nvar of them) */
  uint64_t *srt;  /* their sort keys (bucket, zero L syndrome or not, weight, entry), 2 nvar for the radix sort */
  int nl;         /* the number of rows of L if at most 64 (0: a classical code, or more rows) */
  wi_t width;     /* the words per row of the RW matrix (with nl > 0) */
  word *lrows;    /* with nl > 0: the rows of L (nl rows of `width` words) */
  word *srows;    /* and the rows of L reduced with the pivot rows of a RW step: the L syndromes of the codewords */
} rw_pairs_ws_t;
 
/* The hashes of the columns of L for the RW pairs (see rw_pair_t): the hash of the L syndrome of a vector is the sum of
 * the hashes of its columns.  With at most 64 rows of L, the hash of a column is its mask (the syndrome itself);
 * otherwise, the sum of 64-bit random hashes of its rows (a non-zero syndrome has the hash zero with probability
 * 2^-64).  NULL for a classical code. */
static uint64_t *rw_lhash_init(const params_t * const p) {
  const colmask_t * const cm = p->maskL;
  const int nvar = p->spaH->cols;
  if (p->classical || !cm || cm->cols < nvar) return NULL;
  uint64_t * const h = malloc((size_t)nvar * sizeof(uint64_t));
  if (!h) ERROR("memory allocation");
  if (cm->kw == 1) {
    for (int j = 0; j < nvar; j++) h[j] = cm->mask[j];
  } else {
    uint64_t * const R = malloc((size_t)cm->kw * 64 * sizeof(uint64_t));
    if (!R) ERROR("memory allocation");
    uint64_t st = mix64((uint64_t)p->seed ^ 0x2545f4914f6cdd1dULL);
    for (int t = 0; t < cm->kw * 64; t++) R[t] = splitmix64(&st);
    for (int j = 0; j < nvar; j++) {
      uint64_t s = 0;
      for (int w = 0; w < cm->kw; w++) {
        for (word m = cm->mask[(size_t)j * cm->kw + w]; m; m &= m - 1) s ^= R[w * 64 + __builtin_ctzll(m)];
      }
      h[j] = s;
    }
    free(R);
  }
  return h;
}
 
/* The workspace of the RW pairs of a thread (`ws->ent` NULL without pairs), for the RW matrix `mH` (`nvar` columns) */
static void rw_pairs_ws_init(rw_pairs_ws_t * const ws, const distfork_ctx_t * const ctx, const mzd_t * const mH,
                             const int nvar) {
  memset(ws, 0, sizeof(*ws));
  if (ctx->p->pairs == 0) return;
  ws->ent = malloc((size_t)nvar * sizeof(rw_pair_t));
  ws->srt = malloc((size_t)2 * nvar * sizeof(uint64_t));
  if (!ws->ent || !ws->srt) ERROR("memory allocation");
  const colmask_t * const cm = ctx->p->maskL;
  if (ctx->rw_lhash && cm->kw == 1 && cm->rows > 0) {
    ws->nl = cm->rows;
    ws->width = mH->width;
    ws->lrows = calloc((size_t)ws->nl * ws->width, sizeof(word));
    ws->srows = malloc((size_t)ws->nl * ws->width * sizeof(word));
    if (!ws->lrows || !ws->srows) ERROR("memory allocation");
    for (int j = 0; j < nvar; j++) {
      for (word m = cm->mask[j]; m; m &= m - 1)
        ws->lrows[(size_t)__builtin_ctzll(m) * ws->width + (j >> 6)] |= (word)1 << (j & 63);
    }
  }
}
 
static void rw_pairs_ws_free(rw_pairs_ws_t * const ws) {
  free(ws->ent);
  free(ws->srt);
  free(ws->lrows);
  free(ws->srows);
  memset(ws, 0, sizeof(*ws));
}
 
/* A codeword of weight `cnt` found by RW is a hit (see rw_hits1, rw_hits2) if not heavier than the upper bound */
static inline bool rw_is_hit(distfork_ctx_t * const ctx, const int cnt) {
  const int dmax = atomic_load_explicit(&ctx->dmax, memory_order_relaxed);
  return (dmax <= 0) || (cnt <= dmax);
}
 
/* Weight limit (exclusive) of the codewords of the RW pairs: as for the single columns (rw_cw_limit()), but at most the
 * upper bound plus one (plus dW with outC or maxC), and 0 without an upper bound: the pairs would otherwise give far
 * too many candidates (of any weight, or the heavier codewords of the `min_hits` statistic) */
static inline int rw_pairs_limit(distfork_ctx_t * const ctx) {
  const params_t * const p = ctx->p;
  const int dmax = atomic_load_explicit(&ctx->dmax, memory_order_relaxed);
  if (dmax <= 0) return 0;
  long long int cap = (long long int)dmax + 1;
  if ((p->outC != NULL || p->maxC > 0) && p->dW > 0) cap += p->dW;
  const int lim = rw_cw_limit(ctx);
  return (lim < cap) ? lim : (int)cap;
}
 
/* The pairs visited between the polls of the stop flags and of the upper bound in rw_pairs_step(): with few rows l, a
 * bucket can hold many columns (e.g., with pairs=1, two buckets of (n - r)/2 columns each), so that a step would
 * otherwise not end in time at a timeout */
#define RW_PAIRS_POLL 65536
 
/* The pairs of a RW step (see rw_pair_t), after its single columns: the reduced matrix `mH` (Gauss-Jordan: `rank` pivot
 * rows, whose pivot columns are `pivs`) and its transpose `mHT` (row c: the column c in the pivot rows), `l` pivot rows
 * of the buckets, the workspace `ws`, and the support `ee`; `src` as in rw_record_codeword().  Adds the number of the
 * pairs checked to `*checked`, and of the hits to `*hits` (see rw_is_hit()); returns l (0: no pairs).  The bucket of a
 * column is a hash of its l bits (64 bits; the sort uses its top bits: the buckets with different bits are merged
 * rarely, which only adds pairs). */
static int rw_pairs_step(distfork_ctx_t * const ctx, mzd_t * const mH, mzd_t * const mHT,
                         const rci_t * const pivs, const word * const piv_mask, const int rank, int l,
                         rw_pairs_ws_t * const ws, rci_t * const ee, uint64_t * const rng_state, const int tid,
                         const int src, long long * const checked, long long * const hits) {
  params_t * const p = ctx->p;
  const colmask_t * const maskL = p->maskL;
  const int classical = p->classical;
  const uint64_t * const lhash = ctx->rw_lhash;
  const int nvar = p->spaH->cols;
  rw_pair_t * const ent = ws->ent;
  uint64_t * const srt = ws->srt;
  int lim = rw_pairs_limit(ctx);
  if (lim < 3 || l < 1 || l > rank) return 0; /* (the weight of a pair is at least 2) */
  /* the sort keys: the top kb bits of the bucket's hash, the bit "the L syndrome is not known to be zero", the weight
   * (wb bits), and the entry (ib bits) */
  int wb = 1, ib = 1;
  while ((1L << wb) <= (long)rank) wb++;
  while ((1L << ib) < (long)(nvar - rank)) ib++;
  const int kb = (l + 8 < 63 - wb - ib) ? l + 8 : 63 - wb - ib;
  if (kb < 1) return 0;
  const int sh_grp = ib + wb, sh_key = sh_grp + 1;
  const uint64_t wmask = ((uint64_t)1 << wb) - 1, imask = ((uint64_t)1 << ib) - 1;
  const int aw = (rank + 63) >> 6;
  const word last_mask = (rank & 63) ? (((word)1 << (rank & 63)) - 1) : ~(word)0; /* bits of rows < rank */
  /* l distinct random pivot rows, as the masks `kmask` of the words `kword` of the rows */
  int nk = 0, kword[RW_PAIRS_MAX_BITS];
  word kmask[RW_PAIRS_MAX_BITS];
  for (int i = 0; i < l; i++) {
    int r, j;
    do {
      r = rand_uniform_thread(rank, rng_state);
      for (j = 0; j < nk && kword[j] != (r >> 6); j++) continue;
    } while (j < nk && ((kmask[j] >> (r & 63)) & 1));
    if (j == nk) {
      kword[nk] = r >> 6;
      kmask[nk++] = 0;
    }
    kmask[j] |= (word)1 << (r & 63);
  }
  /* the L syndromes of the codewords of all columns (with at most 64 rows of L): the rows of L reduced with the pivot
   * rows (`lhash` is then the column masks of L) */
  const int nl = ws->nl;
  const wi_t width = ws->width;
  if (nl > 0) {
    memcpy(ws->srows, ws->lrows, (size_t)nl * width * sizeof(word));
    for (int i = 0; i < rank; i++) {
      word g = lhash[pivs[i]];
      if (!g) continue;
      const word * const ri = mzd_row(mH, i);
      for (; g; g &= g - 1) rw_row_xor(ws->srows + (size_t)__builtin_ctzll(g) * width, ri, width);
    }
  }
  const int wsig = lim + 1; /* (otherwise) the L syndromes of the columns up to this weight */
  int K = 0;
  for (int col = 0; col < nvar; col++) {
    if ((piv_mask[col >> 6] >> (col & 63)) & 1) continue;
    const word * const row = mzd_row(mHT, col);
    int wt = 0;
    for (int iw = 0; iw < aw - 1; iw++) wt += __builtin_popcountll(row[iw]);
    wt += __builtin_popcountll(row[aw - 1] & last_mask);
    uint64_t key = 0;
    for (int j = 0; j < nk; j++) key = (key ^ (row[kword[j]] & kmask[j])) * 0x9e3779b97f4a7c15ULL;
    key >>= 64 - kb;
    uint64_t sig = 0;
    int known = 0;
    if (nl > 0) {
      const word * const sc = ws->srows + (col >> 6);
      for (int t = 0; t < nl; t++) sig |= ((sc[(size_t)t * width] >> (col & 63)) & 1) << t;
      known = 1;
    } else if (lhash && wt <= wsig) {
      sig = lhash[col];
      for (int iw = 0; iw < aw; iw++) {
        for (word x = (iw < aw - 1) ? row[iw] : (row[iw] & last_mask); x; x &= x - 1)
          sig ^= lhash[pivs[(iw << 6) + __builtin_ctzll(x)]];
      }
      known = 1;
    }
    const uint64_t grp = !known || (sig != 0);
    ent[K].row = row;
    ent[K].sig = sig;
    ent[K].col = col;
    ent[K].known = known;
    srt[K] = ((uint64_t)key << sh_key) | (grp << sh_grp) | ((uint64_t)wt << ib) | (uint64_t)K;
    K++;
  }
  rw_radix_sort(srt, srt + nvar, K, sh_key + kb);
 
  /* with one row of L (known syndromes), the other columns have the same syndrome: only their pairs with the columns
   * with zero syndromes */
  const bool other_pairs = (nl != 1);
  const bool irred = (p->irred != 0);
  long long int nchk = 0, nhits = 0, nvis = 0, next_poll = RW_PAIRS_POLL; /* nvis: the pairs visited, also skipped */
  long long int nred = 0; /* irred=1: the pairs of disjoint circuits skipped (reducible sums) */
  for (int b0 = 0; b0 < K;) {
    const uint64_t bkey = srt[b0] >> sh_key;
    int b1 = b0 + 1;
    while (b1 < K && (srt[b1] >> sh_key) == bkey) b1++;
    if (b1 - b0 >= 2) {
      if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed) ||
          atomic_load_explicit(&ctx->rw_stop_flag, memory_order_relaxed)) break;
      lim = rw_pairs_limit(ctx); /* (the upper bound may have decreased) */
      if (lim < 3) break;
      int lx = lim - 3; /* the largest weight of the sum of the two columns */
      int g0 = b0;      /* [b0, g0): the columns with zero L syndromes; [g0, b1): the others */
      while (g0 < b1 && !((srt[g0] >> sh_grp) & 1)) g0++;
      int lo = g0;
      const int x_end = other_pairs ? b1 - 1 : g0;
      for (int x = b0; x < x_end; x++) {
        if (nvis >= next_poll) { /* a bucket can be large (few rows l): the stop flags and the upper bound, also here */
          next_poll = nvis + RW_PAIRS_POLL;
          if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed) ||
              atomic_load_explicit(&ctx->rw_stop_flag, memory_order_relaxed)) goto done;
          lim = rw_pairs_limit(ctx);
          if (lim < 3) goto done;
          lx = lim - 3;
        }
        const int wx = (int)((srt[x] >> ib) & wmask);
        const rw_pair_t * const ex = ent + (srt[x] & imask);
        int y = x + 1;
        if (x < g0) { /* (a column with a zero L syndrome: only the pairs with the others) */
          while (lo < b1 && (int)((srt[lo] >> ib) & wmask) < wx - lx) lo++;
          y = lo;
        }
        const int y0 = y;
        for (; y < b1; y++) {
          const int wy = (int)((srt[y] >> ib) & wmask);
          if (wy - wx > lx) break;
          const rw_pair_t * const ey = ent + (srt[y] & imask);
          if (x >= g0 && ex->known && ey->known && ex->sig == ey->sig) continue;
          const word * const ra = ex->row;
          const word * const rb = ey->row;
          nchk++;
          int c = 0;
          for (int iw = 0; iw < aw - 1 && c <= lx; iw++) c += __builtin_popcountll(ra[iw] ^ rb[iw]);
          if (c <= lx) c += __builtin_popcountll((ra[aw - 1] ^ rb[aw - 1]) & last_mask);
          if (c > lx) continue;
          /* irred=1: the sum of the two circuits is reducible iff they are disjoint, i.e., the weight of the sum of
           * their pivot parts is wx + wy (a kernel vector on the pivots and the two columns is a sum of the two
           * circuits, so that the circuits are the only candidates for a proper subset with zero syndrome) */
          if (irred && c == wx + wy) {
            nred++;
            continue;
          }
          int cnt = 0;
          ee[cnt++] = ex->col;
          ee[cnt++] = ey->col;
          for (int iw = 0; iw < aw; iw++) {
            for (word z = (ra[iw] ^ rb[iw]) & ((iw < aw - 1) ? ~(word)0 : last_mask); z; z &= z - 1)
              ee[cnt++] = pivs[(iw << 6) + __builtin_ctzll(z)];
          }
          if (classical || colmask_syndrome_non_zero(maskL, cnt, ee)) {
            rci_quick_sort(ee, cnt);
            if (rw_is_hit(ctx, cnt)) nhits++;
            rw_record_codeword(ctx, ee, cnt, tid, 0, src);
            lim = rw_pairs_limit(ctx);
            if (lim < 3) goto done;
            lx = lim - 3;
          }
        }
        nvis += y - y0;
      }
    }
    b0 = b1;
  }
done:
  *checked += nchk;
  *hits += nhits;
  if (nred > 0) atomic_fetch_add_explicit(&ctx->cw_red_rw, nred, memory_order_relaxed);
  return l;
}
 
/* Run RW batch of the global steps step0 .. step0+n_steps-1; returns the number of completed steps.  `el`: the thread's
 * RW elimination (blocked or column by column, see rw_elim_t); `pws`: the workspace of the RW pairs (see
 * rw_pairs_ws_init()) */
static int run_rw_steps(distfork_ctx_t *ctx, long step0, int n_steps,
                        mzd_t *mH, mzd_t *mHT, rci_t *ee,
                        mzp_t *perm, mzp_t *pivs, word *piv_mask, rw_elim_t *el, rw_pairs_ws_t *pws,
                        int *eff_nrows_ptr,
                        int *visited_cols, int *visited_checks, int *col_queue,
                        int *visit_marker, uint64_t *rng_state, int tid) {
  params_t * const p = ctx->p;
  const colmask_t * const maskL = p->maskL; /* NULL for a classical code */
  const int nvar = p->spaH->cols;
  const int classical = p->classical;
  const int kwin = p->kwin;
  const int win_mode = p->win_mode;
  int eff_nrows = *eff_nrows_ptr;
  int n_done = 0;
  long n_unif = 0; /* completed steps with uniform random permutations (for the information-set estimate) */
  long long int n_pairs = 0, hits1 = 0, hits2 = 0; /* see rw_pairs_checked, rw_hits1, rw_hits2 */
  double t_last = get_time_sec(), acc_t = 0.0;
  int acc_n = 0;
 
  for (int step = 0; step < n_steps; step++) {
    if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed) ||
        atomic_load_explicit(&ctx->rw_stop_flag, memory_order_relaxed)) {
      break;
    }
 
    int eff_kwin = kwin;
    int eff_win_mode = win_mode;
    if (eff_kwin == 0 && nvar >= 500 && ((step0 + step) & 1) == 1) { /* hybrid: every other global step */
      eff_kwin = hybrid_window(p, nvar);
      eff_win_mode = 1;
    }
    const bool uniform = !(eff_kwin > 0 && eff_kwin < nvar);
 
    if (eff_kwin > 0 && eff_kwin < nvar) {
      int seed_col = rand_uniform_thread(nvar, rng_state);
      localized_window_perm(perm, nvar, seed_col, eff_kwin, eff_win_mode,
                            p->spaH, ctx->mHT_cc,
                            visited_cols, visited_checks, col_queue,
                            visit_marker, rng_state);
    } else {
      pivs = mzp_rand_thread(pivs, nvar, rng_state);
      mzp_set_ui(perm, 1);
      perm = perm_p_trans(perm, pivs, 0);
    }
 
    memset(piv_mask, 0, mH->width * sizeof(word));
    /* the elimination (see rw_elim_t): column by column, or blocked */
    int timed = -1; /* tuning: time this step's elimination, column by column (0) or blocked (1) */
    if (el->choice < 0 && el->steps > 0) timed = el->steps & 1;
    const double t_elim = (timed >= 0) ? get_time_sec() : 0.0;
    int rank = 0;
    if (el->choice == 0 || timed == 0) {
      for (int i = 0; i < nvar && rank < eff_nrows; i++) {
        /* a step on a large matrix can take many seconds: check the stop flags also every 16 columns */
        if ((i & 15) == 15 && (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed) ||
                               atomic_load_explicit(&ctx->rw_stop_flag, memory_order_relaxed))) {
          rank = -1;
          break;
        }
        int col = perm->values[i];
        if (gauss_one_rows(mH, col, rank, eff_nrows)) {
          pivs->values[rank++] = col;
          piv_mask[col >> 6] |= (word)1 << (col & 63);
        }
      }
    } else {
      rank = rw_eliminate(ctx, mH, eff_nrows, nvar, perm->values, pivs->values, piv_mask, el);
    }
    /* abandon the step: `eff_nrows` is kept, as rows rank..eff_nrows-1 are not processed (same row space) */
    if (rank < 0) break;
    if (el->choice < 0) {
      if (timed >= 0) el->t[timed] += get_time_sec() - t_elim;
      if (++el->steps > 2 * RW_TUNE_STEPS) {
        el->choice = (el->t[1] < el->t[0]) ? 1 : 0;
        if (tid == 0 && (p->debug & DBG_TIMING)) {
          fprintf(stderr, "# [thread 0] RW elimination: %s (blocked: %.2fx the time of column by column)\n",
                  el->choice ? "blocked" : "column by column", el->t[1] / (el->t[0] > 0.0 ? el->t[0] : 1e-9));
        }
      }
    }
    eff_nrows = rank;
    *eff_nrows_ptr = eff_nrows;
    if (atomic_load_explicit(&ctx->rw_rank, memory_order_relaxed) != rank) { /* rank of H (the same in every step) */
      atomic_store_explicit(&ctx->rw_rank, rank, memory_order_relaxed);
    }
 
    mzd_transpose(mHT, mH);
 
    const int active_width = (rank + 63) >> 6;
    const word last_mask = (rank & 63) ? (((word)1 << (rank & 63)) - 1) : ~(word)0; /* bits of rows < rank */
    for (int col = 0; col < nvar; col++) {
      if ((piv_mask[col >> 6] >> (col & 63)) & 1) continue;
      const int limit = rw_cw_limit(ctx);
      word *rawrow = mzd_row(mHT, col);
      /* the weight 1 + (the number of pivot rows with a 1 in column col) first: word-wise popcount, which ends as soon
       * as the weight reaches the limit (most candidates are much heavier) */
      int wt = 1;
      for (int iw = 0; iw < active_width - 1 && wt < limit; iw++) wt += __builtin_popcountll(rawrow[iw]);
      if (wt < limit && active_width > 0) wt += __builtin_popcountll(rawrow[active_width - 1] & last_mask);
      if (wt >= limit) continue;
 
      int cnt = 0;
      ee[cnt++] = col;
      rci_t j = -1;
      while (cnt < wt) {
        j = nextelement(rawrow, active_width, j);
        if (j == -1 || j >= rank) break;
        ee[cnt++] = pivs->values[j++];
      }
 
      /* the check L c != 0 needs no sorted support: only the codewords recorded are sorted (hash key) */
      if (classical || colmask_syndrome_non_zero(maskL, cnt, ee)) {
        rci_quick_sort(ee, cnt);
        if (rw_is_hit(ctx, cnt)) hits1++;
        rw_record_codeword(ctx, ee, cnt, tid, 0, uniform ? CW_SRC_RW_UNIFORM : CW_SRC_RW_WINDOW);
      }
    }
    if (pws->ent) { /* the pairs of non-pivot columns (p=2) */
      int l = rw_pairs_bits(p, nvar, rank, atomic_load_explicit(&ctx->dmax, memory_order_relaxed));
      if (l > 0)
        l = rw_pairs_step(ctx, mH, mHT, pivs->values, piv_mask, rank, l, pws, ee, rng_state, tid,
                          uniform ? CW_SRC_RW_UNIFORM : CW_SRC_RW_WINDOW, &n_pairs, &hits2);
      if (l > 0 && atomic_load_explicit(&ctx->rw_pairs_l, memory_order_relaxed) != l)
        atomic_store_explicit(&ctx->rw_pairs_l, l, memory_order_relaxed);
    }
    atomic_fetch_add(&ctx->rw_steps_completed, 1);
    n_done++;
    if (uniform) n_unif++;
    rw_time_step(ctx, &t_last, &acc_t, &acc_n, n_done == 1);
  }
  rw_time_publish(ctx, &acc_t, &acc_n);
  atomic_fetch_add_explicit(&ctx->rw_unif_steps, n_unif, memory_order_relaxed);
  if (n_pairs > 0) atomic_fetch_add_explicit(&ctx->rw_pairs_checked, n_pairs, memory_order_relaxed);
  if (hits1 > 0) atomic_fetch_add_explicit(&ctx->rw_hits1, hits1, memory_order_relaxed);
  if (hits2 > 0) atomic_fetch_add_explicit(&ctx->rw_hits2, hits2, memory_order_relaxed);
  return n_done;
}
 
/* After a completed RW step, only the first `rank` = rank(H) rows of the reduced thread-local matrix `*mH` are non-zero
 * (they span the rows of H), and the later steps eliminate only these rows: keep just these rows, so that the
 * transposition in each RW step skips the zero rows (for H with redundant rows) */
static void rw_shrink_matrices(mzd_t **mH, mzd_t **mHT, const int rank) {
  pthread_mutex_lock(&m4ri_mem_mutex);
  mzd_t * const top = mzd_init_window(*mH, 0, 0, rank, (*mH)->ncols);
  mzd_t * const small = mzd_copy(NULL, top);
  mzd_free_window(top);
  mzd_free(*mH);
  *mH = small;
  mzd_free(*mHT);
  *mHT = mzd_init(small->ncols, rank);
  pthread_mutex_unlock(&m4ri_mem_mutex);
}
 
/* Whether row `r` of N has a non-zero entry in a column of the current RW window (visited_cols[j] == marker) */
static inline bool ksub_row_in_window(const mzd_t * const N, const int r, const int nvar,
                                      const int * const visited_cols, const int marker) {
  const word * const raw_n = mzd_row_cons(N, r);
  for (int j = nextelement(raw_n, N->width, -1); j >= 0 && j < nvar; j = nextelement(raw_n, N->width, j + 1)) {
    if (visited_cols[j] == marker) return true;
  }
  return false;
}
 
/* Run compact subspace RW batch (ksub > 0) of the global steps step0 .. step0+n_steps-1; returns the number of
 * completed steps.  `nidx` is the thread's permutation of the row indices 0..nu-1 of N, used to sample ksub distinct
 * rows (partial Fisher-Yates shuffle) */
static int run_rw_steps_ksub(distfork_ctx_t *ctx, long step0, int n_steps,
                             mzd_t *M_sub, rci_t *ee, mzp_t *perm, mzp_t *pivs, int *nidx,
                             int *visited_cols, int *visited_checks, int *col_queue,
                             int *visit_marker, uint64_t *rng_state, int tid) {
  params_t * const p = ctx->p;
  const colmask_t * const maskL = p->maskL; /* NULL for a classical code */
  const int nvar = p->spaH->cols;
  const int classical = p->classical;
  const int nu = ctx->nu;
  const int ksub = M_sub->nrows;
  const int kwin = p->kwin;
  const int win_mode = p->win_mode;
 
  if (nu <= 0 || ksub <= 0 || !ctx->N_global) { /* not reached: no subspace to sample, count the steps as done */
    atomic_fetch_add(&ctx->rw_steps_completed, n_steps);
    return 0;
  }
 
  int n_done = 0;
  double t_last = get_time_sec(), acc_t = 0.0;
  int acc_n = 0;
  pthread_rwlock_rdlock(&ctx->basis_rwlock);
  const mzd_t * const N = ctx->N_global;
 
  for (int step = 0; step < n_steps; step++) {
    if (atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed) ||
        atomic_load_explicit(&ctx->rw_stop_flag, memory_order_relaxed)) {
      break;
    }
 
    /* 1. Generate column permutation (localized window or uniform) */
    int eff_kwin = kwin;
    int eff_win_mode = win_mode;
    if (eff_kwin == 0 && nvar >= 500 && ((step0 + step) & 1) == 1) { /* hybrid: every other global step */
      eff_kwin = hybrid_window(p, nvar);
      eff_win_mode = 1;
    }
 
    if (eff_kwin > 0 && eff_kwin < nvar) {
      int seed_col = rand_uniform_thread(nvar, rng_state);
      localized_window_perm(perm, nvar, seed_col, eff_kwin, eff_win_mode,
                            p->spaH, ctx->mHT_cc,
                            visited_cols, visited_checks, col_queue,
                            visit_marker, rng_state);
    } else {
      pivs = mzp_rand_thread(pivs, nvar, rng_state);
      mzp_set_ui(perm, 1);
      perm = perm_p_trans(perm, pivs, 0);
    }
 
    /* 2. Sample ksub distinct rows of N into M_sub: partial Fisher-Yates shuffle of the row indices (rows drawn
     * with replacement would repeat and reduce the rank), with up to 4 redraws to prefer rows overlapping the window */
    const int marker = *visit_marker;
    const bool in_window = (eff_kwin > 0 && eff_kwin < nvar && visited_cols);
    for (int i = 0; i < ksub; i++) {
      int j = i + rand_uniform_thread(nu - i, rng_state);
      for (int attempt = 0; in_window && attempt < 4 && !ksub_row_in_window(N, nidx[j], nvar, visited_cols, marker);
           attempt++) {
        j = i + rand_uniform_thread(nu - i, rng_state);
      }
      SWAPINT(nidx[i], nidx[j]);
      mzd_copy_row(M_sub, i, N, nidx[i]);
    }
 
    /* 3. Echelonize M_sub directly in L1/L2 cache using permuted column order */
    int rank = 0;
    for (int i = 0; i < nvar && rank < ksub; i++) {
      int col = perm->values[i];
      if (gauss_one(M_sub, col, rank)) {
        rank++;
      }
    }
 
    /* 4. Extract candidate codewords directly from reduced rows (already in sorted order) */
    for (int ir = 0; ir < rank; ir++) {
      int cnt = 0;
      const int limit = rw_cw_limit(ctx);
 
      word *rawrow = mzd_row(M_sub, ir);
      rci_t j = -1;
      while (cnt < limit) {
        j = nextelement(rawrow, M_sub->width, j);
        if (j == -1 || j >= nvar) break;
        ee[cnt++] = j++;
      }
 
      if (cnt > 0 && cnt < limit && (classical || colmask_syndrome_non_zero(maskL, cnt, ee))) {
        /* irred=1: a reduced row of the sampled subspace is not a circuit, so it can be reducible */
        if (!cw_irreducible(ctx, ctx->wargs[tid].irr_ws, cnt, ee)) {
          atomic_fetch_add_explicit(&ctx->cw_red_rw, 1, memory_order_relaxed);
          continue;
        }
        rw_record_codeword(ctx, ee, cnt, tid, ksub,
                           (eff_kwin > 0 && eff_kwin < nvar) ? CW_SRC_RW_WINDOW : CW_SRC_RW_UNIFORM);
      }
    }
    atomic_fetch_add(&ctx->rw_steps_completed, 1);
    n_done++;
    rw_time_step(ctx, &t_last, &acc_t, &acc_n, n_done == 1);
  }
  rw_time_publish(ctx, &acc_t, &acc_n);
 
  pthread_rwlock_unlock(&ctx->basis_rwlock);
 
  /* 5. Periodic adaptive basis refresh if enabled */
  if (p->refresh > 0 && !atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed) &&
      !atomic_load_explicit(&ctx->rw_stop_flag, memory_order_relaxed)) {
    long done = atomic_load_explicit(&ctx->rw_steps_completed, memory_order_relaxed);
    long next_ref = atomic_load_explicit(&ctx->next_refresh_step, memory_order_relaxed);
    if (done >= next_ref && next_ref > 0) {
      if (atomic_compare_exchange_strong(&ctx->next_refresh_step, &next_ref,
                                         done + p->refresh)) {
        int cw_wt = 0;
        pthread_mutex_lock(&ctx->cw_mutex);
        if (p->codewords && p->min_w < nvar) {
          cw_vec_t *cw, *tmp;
          HASH_ITER(hh, p->codewords, cw, tmp) {
            if (cw->weight == p->min_w) {
              cw_wt = cw->weight;
              for (int k = 0; k < cw_wt; k++) ee[k] = cw->arr[k];
              break;
            }
          }
        }
        pthread_mutex_unlock(&ctx->cw_mutex);
 
        pthread_rwlock_wrlock(&ctx->basis_rwlock);
        pthread_mutex_lock(&m4ri_mem_mutex);
        refresh_nullspace_basis(ctx->N_global, cw_wt > 0 ? ee : NULL, cw_wt, rng_state);
        pthread_mutex_unlock(&m4ri_mem_mutex);
        pthread_rwlock_unlock(&ctx->basis_rwlock);
      }
    }
  }
  return n_done;
}
 
/* Worker thread main loop */
static void *worker_thread_func(void *arg) {
  worker_arg_t *warg = (worker_arg_t *)arg;
  distfork_ctx_t *ctx = warg->ctx;
  int tid = warg->tid;
  const int nvar = ctx->p->spaH->cols;
  /* only workers 0..rw_threads-1 run RW (method 1 or 3) and allocate dense RW matrices; all workers may run CC */
  const bool enable_rw = ((ctx->p->method & 1) != 0) && (tid < ctx->rw_threads);
  const bool use_ksub = enable_rw && (ctx->p->ksub > 0);
  const int ksub_eff = (use_ksub && ctx->nu > 0) ? minint(ctx->p->ksub, ctx->nu) : 0;
  /* (warg->min_swei is initialized by the main thread before the workers are launched) */
 
  /* Thread-local RW matrices (allocated safely only if RW is enabled) */
  mzd_t *mH = NULL;
  mzd_t *mHT_rw = NULL;
  mzd_t *M_sub = NULL;
  int *nidx = NULL; /* ksub: permutation of the row indices of N, for sampling distinct rows */
  rci_t *ee = NULL;
  mzp_t *perm = NULL;
  mzp_t *pivs = NULL;
  word *piv_mask = NULL;
  rw_elim_t rw_el = {NULL, 0, 0, 0, 0, {0.0, 0.0}, 0}; /* the RW elimination of this thread */
  rw_pairs_ws_t rw_pws; /* the RW pairs of this thread (see rw_pairs_ws_init()) */
  memset(&rw_pws, 0, sizeof(rw_pws));
  int eff_nrows = ctx->p->spaH->rows;
  int *visited_cols = NULL;
  int *visited_checks = NULL;
  int *col_queue = NULL;
  int visit_marker = 0;
  /* RNG stream of this thread: a hashed starting point of splitmix64 (consecutive starting points would give
   * shifted copies of one stream) */
  uint64_t rng_state = mix64((uint64_t)ctx->p->seed ^ mix64((uint64_t)tid + 0x517cc1b727220a95ULL));
 
  if (enable_rw) {
    ee = malloc((nvar + 2) * sizeof(rci_t));
    perm = safe_mzp_init(nvar);
    pivs = safe_mzp_init(nvar);
    if (ctx->p->kwin > 0 || nvar >= 500) {
      visited_cols = calloc(nvar, sizeof(int));
      visited_checks = calloc(ctx->p->spaH->rows, sizeof(int));
      col_queue = calloc(nvar, sizeof(int));
    }
    if (use_ksub) {
      if (ksub_eff > 0) {
        M_sub = safe_mzd_init(ksub_eff, nvar);
        nidx = malloc((size_t)ctx->nu * sizeof(int));
        if (!nidx) ERROR("memory allocation");
        for (int i = 0; i < ctx->nu; i++) nidx[i] = i;
      }
    } else {
      mH = safe_mzd_from_csr(NULL, ctx->p->spaH);
      mHT_rw = safe_mzd_init(nvar, ctx->p->spaH->rows);
      piv_mask = calloc(mH->width, sizeof(word));
      /* the tables of rw_eliminate(), with its parameters for this matrix (64-byte aligned rows) */
      rw_el.tstride = ((size_t)mH->width + 7) & ~(size_t)7;
      rw_block_params(mH->nrows, mH->width, &rw_el.k, &rw_el.nt);
      rw_el.choice = (rw_el.nt == 0) ? 0 : ((mH->nrows >= 512) ? 1 : -1); /* (tuning between 128 and 512 rows) */
      const size_t tbytes = (rw_el.nt > 0 ? ((size_t)rw_el.nt << rw_el.k) : 1) * rw_el.tstride * sizeof(word);
      rw_el.table = aligned_alloc(64, tbytes);
      if (!piv_mask || !rw_el.table) ERROR("memory allocation");
      memset(rw_el.table, 0, tbytes);
      rw_pairs_ws_init(&rw_pws, ctx, mH, nvar);
    }
  }
 
  /* Thread-local CC memory (only for CC, method 2 or 3) */
  const int wmax_alloc = MAX_W - 1;
  one_vec_t *err = NULL, *urr = NULL, **syn = NULL;
  cc_walk_t walk; /* the CC walk with bitset syndromes (smax = 0: ctx->cc_hash, see cc_walk()) */
  memset(&walk, 0, sizeof(walk));
  if (ctx->p->method >= 2) {
    err = calloc(1, sizeof(one_vec_t) + sizeof(int) * (wmax_alloc + 2));
    urr = calloc(1, sizeof(one_vec_t) + sizeof(int) * (wmax_alloc + 2));
    syn = calloc(wmax_alloc + 3, sizeof(one_vec_t *));
    if (!err || !urr || !syn) ERROR("memory allocation");
    for (int i = 0; i <= wmax_alloc + 2; i++) {
      syn[i] = calloc(1, sizeof(one_vec_t) + sizeof(int) * (ctx->p->spaH->rows + 1));
      if (!syn[i]) ERROR("memory allocation");
    }
    if (ctx->cc_hash) {
      walk.ctx = ctx;
      walk.warg = warg;
      walk.mH = ctx->p->spaH;
      walk.mHT = ctx->mHT_cc;
      walk.hs = ctx->cc_hash;
      walk.m = ctx->p->spaH->rows;
      walk.max_col_wt = ctx->max_col_W;
      walk.S = calloc((size_t)(walk.m / 64 + 2), sizeof(uint64_t));
      walk.E = calloc((size_t)(nvar / 64 + 2), sizeof(uint64_t));
      walk.cols = calloc(wmax_alloc + 3, sizeof(int));
      walk.cw = calloc(1, sizeof(one_vec_t) + sizeof(int) * (wmax_alloc + 3));
      if (!walk.S || !walk.E || !walk.cols || !walk.cw) ERROR("memory allocation");
      walk.prune = ctx->p->cc_prune;
      if (walk.prune > 0) { /* the branch and bound, see cc_bnb() */
        const int B = (ctx->cc_B > 0) ? ctx->cc_B : 1;
        walk.max_row_wt = B;
        walk.acnt_col = INT_MIN;
        walk.X = calloc((size_t)(nvar / 64 + 2), sizeof(uint64_t));
        walk.acnt = calloc((size_t)walk.m + 1, sizeof(int));
        walk.kfirst = calloc((size_t)walk.m + 1, sizeof(int));
        walk.fstk = malloc(((size_t)nvar + 1) * sizeof(int));
        walk.brl = malloc((size_t)(wmax_alloc + 3) * (size_t)B * sizeof(int));
        walk.ulist = malloc(((size_t)walk.m + 1) * sizeof(int));
        walk.uacnt = malloc(((size_t)walk.m + 1) * sizeof(int));
        walk.ord = malloc(((size_t)walk.m + 1) * sizeof(int));
        walk.bcnt = malloc(((size_t)B + 2) * sizeof(int));
        walk.mark = calloc((size_t)nvar + 1, sizeof(uint32_t));
        if (!walk.X || !walk.acnt || !walk.kfirst || !walk.fstk || !walk.brl || !walk.ulist || !walk.uacnt ||
            !walk.ord || !walk.bcnt || !walk.mark)
          ERROR("memory allocation");
        if (walk.prune >= 2) {
          walk.cv = malloc(((size_t)nvar + 1) * sizeof(int));
          walk.touched = malloc(((size_t)nvar + 1) * sizeof(int));
          walk.hist = malloc(((size_t)walk.max_col_wt + 2) * sizeof(int));
          if (!walk.cv || !walk.touched || !walk.hist) ERROR("memory allocation");
        }
      }
    }
  }
 
  /* irred=1: the workspace of the irreducibility test of the codewords found by CC and by the experimental ksub RW */
  if (ctx->p->irred && (ctx->p->method >= 2 || use_ksub)) warg->irr_ws = cols_ws_init(ctx->p->spaH->rows);
 
  while (!atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed)) {
    if (ctx->timeout > 0.0 && (get_time_sec() - ctx->start_time >= ctx->timeout)) {
      ctx_stop(ctx, STOP_TIMEOUT);
      break;
    }
 
    bool did_work = false;
 
    /* 1. Try to take CC work if CC is active (method 2 or 3) */
    if (ctx->p->method >= 2 && atomic_load_explicit(&ctx->cc_round_active, memory_order_acquire)) {
      int active = atomic_load_explicit(&ctx->cc_active_workers, memory_order_relaxed);
      int target = atomic_load_explicit(&ctx->cc_target_workers, memory_order_relaxed);
      /* A worker without RW work joins any CC round: method 2, a CC-only worker (tid >= rw_threads), RW stopped by
       * `min_hits`, or all RW steps already claimed by other workers (it would otherwise idle until the round ends) */
      if (!enable_rw || atomic_load_explicit(&ctx->rw_stop_flag, memory_order_relaxed) ||
          atomic_load_explicit(&ctx->rw_steps_started, memory_order_relaxed) >= ctx->total_rw_steps) {
        target = ctx->num_threads;
      }
      const long long next = atomic_load_explicit(&ctx->cc_next, memory_order_relaxed);
      if (active < target && cc_unpack_idx(next) <= cc_unit_last(ctx, cc_unpack_w(next), cc_unpack_lev(next))) {
        if (atomic_compare_exchange_weak(&ctx->cc_active_workers, &active, active + 1)) {
          const double t_join = get_time_sec();
          while (!atomic_load_explicit(&ctx->stop_flag, memory_order_relaxed)) {
            if (ctx->timeout > 0.0 && (get_time_sec() - ctx->start_time >= ctx->timeout)) {
              ctx_stop(ctx, STOP_TIMEOUT);
              break;
            }
            /* claim the next work unit together with the weight and the fixed levels of its round (see cc_pack()) */
            const long long claim = atomic_fetch_add_explicit(&ctx->cc_next, 1, memory_order_relaxed);
            const int w = cc_unpack_w(claim);
            const int lev = cc_unpack_lev(claim);
            const int unit = cc_unpack_idx(claim);
            if (unit > cc_unit_last(ctx, w, lev)) break;
            /* the start index of the unit, and its branches at the fixed levels: the digits of the unit in radix B */
            const int idx = unit / ctx->cc_R[lev];
            int fix[CC_SPLIT_MAX_LEVELS];
            for (int l = 0, k = unit - idx * ctx->cc_R[lev]; l < lev; l++, k /= ctx->cc_B) fix[l] = k % ctx->cc_B;
            /* with the expert `start` list, `idx` enumerates the listed start columns */
            const int col = (ctx->p->start_num > 0) ? ctx->p->start_list[idx] : idx;
 
            if (walk.S) {
              const int swei = ctx->mHT_cc->p[col + 1] - ctx->mHT_cc->p[col];
              if (w > 1) {
                if (swei > 0 && swei <= (w - 1) * ctx->max_col_W) {
                  if (walk.prune > 0) cc_bnb_unit(&walk, col, idx, w, fix, lev);
                  else cc_walk_unit(&walk, col, w, fix, lev);
                }
                continue;
              }
              if (swei > 0) continue;
            }
 
            err->vec[0] = urr->vec[0] = col;
            err->wei = urr->wei = 1;
            syn[1]->wei = 0;
            int swei = walk.S ? 0 : one_csr_row_combine(syn[1], syn[0], ctx->mHT_cc, col);
 
            if (ctx->p->smax && swei > 0 && swei <= ctx->p->smax) {
              if (swei < warg->min_swei[1]) {
                warg->min_swei[1] = swei;
              }
            }
 
            if (w > 1) {
              if (swei > 0 && swei <= (w - 1) * ctx->max_col_W + ctx->p->smax) {
                if (walk.S) {
                  if (walk.prune > 0) cc_bnb_unit(&walk, col, idx, w, fix, lev);
                  else cc_walk_unit(&walk, col, w, fix, lev);
                } else {
                  start_CC_recurs_mt(err, urr, syn, w, ctx->max_col_W,
                                     ctx->p->spaH, ctx->mHT_cc, warg, fix, lev);
                }
              }
            } else {
              if (!swei) {
                int nz = (!ctx->p->maskL) || colmask_syndrome_non_zero(ctx->p->maskL, 1, err->vec);
                if (nz) {
                  params_t * const p = ctx->p;
                  pthread_mutex_lock(&ctx->cw_mutex);
                  p->codewords = codeword_add_maybe(p, err->vec, 1, CW_SRC_OTHER);
                  ctx_update_feed_limit(ctx);
                  if ((p->debug & DBG_CODEWORDS) && atomic_load(&ctx->dmax) != 1) {
                    char prefix[48];
                    snprintf(prefix, sizeof(prefix), "# [thread %d] CC codeword", tid);
                    print_codeword_support(stderr, prefix, err->vec, 1, p->debug);
                  }
                  atomic_store(&ctx->cc_found_weight, 1);
                  atomic_store(&ctx->cc_exact, true); /* weight 1 is always the exact distance */
                  atomic_store(&ctx->dmin, 1);
                  atomic_store(&ctx->dmax, 1);
                  /* d=1 is known; with outC or maxC, the round goes on to collect all codewords of weight 1 */
                  if (codeword_maxc_reached(p)) {
                    ctx_stop(ctx, STOP_MAXC);
                  } else if (!p->outC && p->maxC == 0) {
                    ctx_stop(ctx, STOP_CC_FOUND);
                  }
                  pthread_mutex_unlock(&ctx->cw_mutex);
                }
              }
            }
            err->wei = urr->wei = 0;
          }
          /* CC thread time (added before leaving the round, so the coordinator sees it once the round is over) */
          atomic_fetch_add_explicit(&ctx->cc_busy_ns, (long long)((get_time_sec() - t_join) * 1e9),
                                    memory_order_relaxed);
          atomic_fetch_add_explicit(&ctx->cc_nodes, walk.nodes, memory_order_relaxed);
          walk.nodes = 0;
          atomic_fetch_sub_explicit(&ctx->cc_active_workers, 1, memory_order_release);
          did_work = true;
          continue;
        }
      }
    }
 
    /* 2. Try to take RW work if RW is active (method 1 or 3) */
    if (enable_rw && !atomic_load(&ctx->stop_flag) && !atomic_load(&ctx->rw_stop_flag)) {
      long cur_s = atomic_load(&ctx->rw_steps_started);
      if (cur_s < ctx->total_rw_steps) {
        long batch_size;
        if (ctx->p->chunk_size > 0) {
          batch_size = ctx->p->chunk_size;
        } else {
          int nvar = ctx->p->nvar;
          if (nvar < 500) {
            if (ctx->total_rw_steps >= 50000) batch_size = 500;
            else if (ctx->total_rw_steps >= 1000) batch_size = 250;
            else batch_size = 50;
          } else if (nvar < 5000) {
            if (ctx->total_rw_steps >= 10000) batch_size = 100;
            else batch_size = 50;
          } else {
            /* Large matrices (e.g. n >= 5000): keep chunk bounded to ~2-3s */
            if (ctx->total_rw_steps >= 10000) batch_size = 50;
            else batch_size = 25;
          }
        }
        /* Guided scheduling: a batch takes at most a quarter of the remaining steps per RW thread, so that no thread
         * monopolizes the steps and the RW threads finish together (no idle threads at the end of the steps) */
        if (ctx->rw_threads > 1) {
          const long guided = (ctx->total_rw_steps - cur_s + 4L * ctx->rw_threads - 1) / (4L * ctx->rw_threads);
          if (batch_size > guided) batch_size = (guided > 0) ? guided : 1;
        }
        long target_s = cur_s + batch_size;
        if (target_s > ctx->total_rw_steps) target_s = ctx->total_rw_steps;
        if (atomic_compare_exchange_weak(&ctx->rw_steps_started, &cur_s, target_s)) {
          int n_steps = (int)(target_s - cur_s);
          if (use_ksub) {
            if (ksub_eff > 0) {
              run_rw_steps_ksub(ctx, cur_s, n_steps, M_sub, ee, perm, pivs, nidx,
                                visited_cols, visited_checks, col_queue,
                                &visit_marker, &rng_state, tid);
            } else {
              atomic_fetch_add(&ctx->rw_steps_completed, n_steps);
            }
          } else {
            run_rw_steps(ctx, cur_s, n_steps, mH, mHT_rw, ee, perm, pivs,
                         piv_mask, &rw_el, &rw_pws, &eff_nrows, visited_cols, visited_checks,
                         col_queue, &visit_marker, &rng_state, tid);
            if (eff_nrows > 0 && eff_nrows < mH->nrows) rw_shrink_matrices(&mH, &mHT_rw, eff_nrows);
          }
          did_work = true;
          continue;
        }
      }
    }
 
    if (!did_work) {
      usleep(100);
    }
  }
 
  if (enable_rw) {
    free(piv_mask);
    free(rw_el.table);
    rw_pairs_ws_free(&rw_pws);
    safe_mzp_free(perm);
    safe_mzp_free(pivs);
    free(ee);
    free(visited_cols);
    free(visited_checks);
    free(col_queue);
    free(nidx);
    safe_mzd_free(M_sub);
    safe_mzd_free(mHT_rw);
    safe_mzd_free(mH);
  }
 
  if (syn) {
    for (int i = 0; i <= wmax_alloc + 2; i++) free(syn[i]);
    free(syn);
  }
  free(err);
  free(urr);
  free(walk.S);
  free(walk.E);
  free(walk.cols);
  free(walk.cw);
  free(walk.X);
  free(walk.acnt);
  free(walk.kfirst);
  free(walk.fstk);
  free(walk.brl);
  free(walk.ulist);
  free(walk.uacnt);
  free(walk.ord);
  free(walk.bcnt);
  free(walk.mark);
  free(walk.cv);
  free(walk.touched);
  free(walk.hist);
  warg->irr_ws = cols_ws_free(warg->irr_ws);
 
  return NULL;
}
 
/* Percentage of the start indices of the current (or last) CC round claimed by the workers (in work units) */
static double cc_round_claimed(distfork_ctx_t * const ctx) {
  const long long beg = (long long)ctx->round_beg * ctx->round_R;
  const long long end = ((long long)ctx->round_end + 1) * ctx->round_R; /* one past the last work unit */
  if (end <= beg) return 100.0;
  long long idx = cc_unpack_idx(atomic_load(&ctx->cc_next));
  if (idx > end) idx = end;
  if (idx < beg) idx = beg;
  return 100.0 * (double)(idx - beg) / (double)(end - beg);
}
 
/* State of RW once it no longer runs, for messages */
static const char *rw_end_state(distfork_ctx_t * const ctx) {
  if (atomic_load(&ctx->rw_stop_flag)) return "converged";
  if (atomic_load(&ctx->rw_steps_completed) >= ctx->total_rw_steps) return "all steps done";
  return "stopped";
}
 
/* DBG_STATUS: periodic status line, after 1, 2, 4, ... s, then every 60 s (called by the coordinator only) */
static void ctx_status(distfork_ctx_t * const ctx) {
  params_t * const p = ctx->p;
  if (!(p->debug & DBG_STATUS)) return;
  const double t = get_time_sec() - ctx->start_time;
  if (t < ctx->status_next) return;
  ctx->status_next += ctx->status_dt;
  if (ctx->status_next <= t) ctx->status_next = t + ctx->status_dt;
  ctx->status_dt = (2.0 * ctx->status_dt < 60.0) ? 2.0 * ctx->status_dt : 60.0;
 
  char cc_str[160] = "", rw_str[224] = "";
  if (p->method >= 2) {
    if (ctx->round_open) {
      int len = snprintf(cc_str, sizeof(cc_str), "; CC round w=%d: %.1f%% of start columns claimed, %.3gs",
                         ctx->round_w, cc_round_claimed(ctx), get_time_sec() - ctx->round_t0);
      if (ctx->round_est > 0.0 && ctx->round_threads > 0 && len > 0 && len < (int)sizeof(cc_str))
        snprintf(cc_str + len, sizeof(cc_str) - len, " (predicted %.3gs on %d threads)",
                 ctx->round_est / ctx->round_threads, ctx->round_threads);
    } else {
      snprintf(cc_str, sizeof(cc_str), "; CC %s", ctx->cc_idle ? ctx->cc_idle : "between rounds");
    }
  }
  if ((p->method & 1) && ctx->rw_threads > 0) {
    const long steps = atomic_load(&ctx->rw_steps_completed);
    const double rate = (t > ctx->status_t) ? (double)(steps - ctx->status_steps) / (t - ctx->status_t) : 0.0;
    ctx->status_steps = steps;
    ctx->status_t = t;
    const bool rw_on = !atomic_load(&ctx->rw_stop_flag) && steps < ctx->total_rw_steps;
    int len = snprintf(rw_str, sizeof(rw_str), "; RW: %ld of %ld steps (%.3g/s)%s%s", steps, ctx->total_rw_steps,
                       rate, rw_on ? "" : ", ", rw_on ? "" : rw_end_state(ctx));
    if (p->min_hits > 0 && len > 0 && len < (int)sizeof(rw_str)) {
      cw_hit_stats_t st;
      pthread_mutex_lock(&ctx->cw_mutex);
      codeword_hit_stats(p, &st);
      pthread_mutex_unlock(&ctx->cw_mutex);
      if (st.set_cws > 0)
        snprintf(rw_str + len, sizeof(rw_str) - len, ", <n>=%.2f over %lld cws of weight %d..%d", st.avg,
                 st.set_cws, st.w_lo, st.w_hi);
    }
  }
  fprintf(stderr, "# status %.1fs: bounds [%d, %d]%s%s\n", t, atomic_load(&ctx->dmin), atomic_load(&ctx->dmax),
          cc_str, rw_str);
}
 
/* Start the workers ctx->spawned .. n-1 (main thread only).  If a thread cannot be created, the run goes on with the
 * workers already started, and no more are started. */
static void ctx_spawn(distfork_ctx_t * const ctx, const int n) {
  for (int i = ctx->spawned; i < n && !ctx->spawn_failed; i++) {
    const int err = pthread_create(&ctx->threads[i], NULL, worker_thread_func, &ctx->wargs[i]);
    if (err != 0) {
      if (i == 0) ERROR("cannot create a worker thread (error %d)", err);
      fprintf(stderr, "# Warning: cannot create worker thread %d (error %d), running with %d threads\n", i, err, i);
      ctx->spawn_failed = true;
      break;
    }
    ctx->spawned = i + 1;
  }
}
 
/* Method 2 (unless nothrottle=1) starts only CC_SPAWN_INITIAL workers, and the others for a CC round predicted to take
 * more than CC_SPAWN_WORK thread-seconds per started worker (or of unknown work), or once a round has run for
 * CC_SPAWN_TIME seconds: tiny runs avoid the start-up of many threads, longer runs use all of them (any worker can
 * join a CC round at any time) */
#define CC_SPAWN_INITIAL 4
#define CC_SPAWN_WORK 2e-3
#define CC_SPAWN_TIME 5e-3
 
/* Method 2: some workers not started yet */
static inline bool m2_can_spawn(const distfork_ctx_t * const ctx) {
  return ctx->spawned < ctx->num_threads && !ctx->spawn_failed;
}
 
/* Method 2: start the remaining workers during (or for) the CC round at weight w; `why`: for the message */
static void m2_spawn_rest(distfork_ctx_t * const ctx, const int w, const char * const why) {
  const int n0 = ctx->spawned;
  ctx_spawn(ctx, ctx->num_threads);
  ctx->round_threads = ctx->spawned;
  if ((ctx->p->debug & DBG_TIMING) && ctx->spawned > n0) {
    fprintf(stderr, "# CC round w=%d: started %d more threads, %d in total (%s)\n", w, ctx->spawned - n0, ctx->spawned,
            why);
  }
}
 
/* Method 1 coordinator */
static void run_method1_coordinator(distfork_ctx_t *ctx) {
  while (!atomic_load(&ctx->stop_flag)) {
    if (ctx->timeout > 0.0 && (get_time_sec() - ctx->start_time >= ctx->timeout)) {
      ctx_stop(ctx, STOP_TIMEOUT);
      break;
    }
    if (atomic_load(&ctx->rw_steps_completed) >= ctx->total_rw_steps) {
      ctx_stop(ctx, STOP_RW_STEPS);
      break;
    }
    ctx_status(ctx);
    usleep(1000);
  }
}
 
/* Estimated work of the CC round at weight w in thread-seconds (the unit of `cc_time_per_weight`, i.e., wall time
 * times the number of CC threads): the measured work W(w-1) of the previous round times the growth factor
 * W(w-1)/W(w-2), clamped to [2, 10] (4 if W(w-2) is not known), also returned in `*growth`.  Returns 0 if W(w-1)
 * has not been measured (the first round at a supplied dmin > 1). */
static double cc_work_estimate(const distfork_ctx_t * const ctx, const int w, double * const growth) {
  if (ctx->even_weight_only && !(w & 1)) {
    *growth = 16.0;
    if (w <= 2) return 1e-4;
    if (w > MAX_W) return 0.0;
    const double prev = ctx->cc_time_per_weight[w - 2];
    if (!(prev > 0.0)) return 0.0;
    if (w >= 6 && ctx->cc_time_per_weight[w - 4] > 1e-4) {
      const double g = prev / ctx->cc_time_per_weight[w - 4];
      *growth = (g < 4.0) ? 4.0 : ((g > 100.0) ? 100.0 : g);
    }
    return ((prev > 1e-4) ? prev : 1e-4) * (*growth);
  }
  *growth = 4.0;
  if (w <= 1) return 1e-4; /* one syndrome per column */
  if (w > MAX_W) return 0.0;
  const double prev = ctx->cc_time_per_weight[w - 1];
  if (!(prev > 0.0)) return 0.0;
  if (w >= 3 && ctx->cc_time_per_weight[w - 2] > 1e-4) {
    const double g = prev / ctx->cc_time_per_weight[w - 2];
    *growth = (g < 2.0) ? 2.0 : ((g > 10.0) ? 10.0 : g);
  }
  return ((prev > 1e-4) ? prev : 1e-4) * (*growth);
}
 
/* Number of collected codewords (the collection window, as exported with outC), for messages */
static long long ctx_num_cws(distfork_ctx_t * const ctx) {
  pthread_mutex_lock(&ctx->cw_mutex);
  const long long num = codeword_window_count(ctx->p);
  pthread_mutex_unlock(&ctx->cw_mutex);
  return num;
}
 
/* Method 2 coordinator */
static void run_method2_coordinator(distfork_ctx_t *ctx) {
  const int wmax = ctx->p->wmax;
  const int w_start = ctx->p->noscan ? wmax : (ctx->p->dmin > 1 ? ctx->p->dmin : 1);
  /* with outC or maxC, CC does not stop at the first codeword but collects all codewords of a round */
  const bool collecting = (ctx->p->outC != NULL) || (ctx->p->maxC > 0);
  /* with outC and dW > 0, extra CC rounds export the codewords of weight up to d + dW */
  const int extra_w = (ctx->p->outC && ctx->p->dW > 0) ? ctx->p->dW : 0;
 
  int w_limit = wmax;
  for (int w = w_start; ; w++) {
    if (atomic_load(&ctx->stop_flag)) break;
    /* The current upper bound dmax (given, read from finC, or found by CC) limits the rounds: since the CC kernel
     * caps the cluster weight at dmax + dW, a round w > dmax + extra_w would only enumerate the same clusters
     * again.  Unless codewords are collected, the round w = dmax is not needed once dmin = dmax. */
    const int cur_dmax = atomic_load(&ctx->dmax);
    if (cur_dmax > 0) {
      if (!collecting && atomic_load(&ctx->dmin) >= cur_dmax) {
        if (ctx->p->debug & DBG_PROGRESS) {
          fprintf(stderr, "# bounds coincide: dmin = dmax = %d (CC round w=%d not needed)\n", cur_dmax, w);
        }
        ctx_stop(ctx, STOP_BOUNDS);
        break;
      }
      w_limit = minint(w_limit, cur_dmax + extra_w);
    }
    /* The stop target dstop (if below dmax): no rounds w >= dstop are needed, except the rounds w = dstop..dstop+dW
     * which collect codewords */
    const int dstop = ctx->p->dstop;
    if (dstop > 0 && (cur_dmax == 0 || dstop < cur_dmax)) {
      if (!collecting && atomic_load(&ctx->dmin) >= dstop) {
        if (ctx->p->debug & DBG_PROGRESS) {
          fprintf(stderr, "# lower bound dmin=%d reached dstop=%d (CC round w=%d not needed)\n",
                  atomic_load(&ctx->dmin), dstop, w);
        }
        ctx_stop(ctx, STOP_DSTOP);
        break;
      }
      w_limit = minint(w_limit, collecting ? dstop + extra_w : dstop - 1);
    }
    if (w > w_limit) {
      ctx->stop_w = w_limit;
      ctx_stop(ctx, (dstop > 0 && atomic_load(&ctx->dmin) >= dstop) ? STOP_DSTOP : STOP_CC_DONE);
      break;
    }
    if (ctx->even_weight_only && (w & 1)) {
      if (w == atomic_load(&ctx->dmin)) atomic_store(&ctx->dmin, w + 1);
      continue;
    }
 
    double now = get_time_sec();
    double remaining_time = ctx->timeout - (now - ctx->start_time);
    if (ctx->timeout > 0.0 && remaining_time <= 0.0) {
      ctx_stop(ctx, STOP_TIMEOUT);
      break;
    }
 
    /* With a timeout, a round predicted not to finish in time is started anyway: it can still end early with a
     * codeword of weight w, which gives the exact distance (all lower weights have been analyzed) */
    double growth;
    const double t_cc_work = cc_work_estimate(ctx, w, &growth); /* thread-seconds (0: unknown) */
    if (ctx->timeout > 0.0 && (ctx->p->debug & DBG_PROGRESS)) {
      const double t_cc_est = t_cc_work / ctx->num_threads; /* wall time */
      if (t_cc_est > remaining_time) {
        fprintf(stderr, "# CC for w=%d (est %.2fs) may not finish in the remaining time %.2fs, searching anyway "
                "(dmin=%d)\n", w, t_cc_est, remaining_time, atomic_load(&ctx->dmin));
      }
    }
 
    const int beg = cc_round_beg(ctx->p);
    const int end = cc_round_end(ctx->p, w);
    /* all weights below w analyzed (or a supplied dmin): a codeword of weight w found is exact */
    const bool certified = (w == atomic_load(&ctx->dmin));
    /* the workers not started yet, for a longer round (or one of unknown work) */
    if (m2_can_spawn(ctx) && !(t_cc_work > 0.0 && t_cc_work <= CC_SPAWN_WORK * ctx->spawned)) {
      char why[64];
      if (t_cc_work > 0.0) snprintf(why, sizeof(why), "predicted %.3g thread-s", t_cc_work);
      else snprintf(why, sizeof(why), "work not known");
      m2_spawn_rest(ctx, w, why);
    }
    const int lev = cc_round_levels(ctx, w, beg, end, t_cc_work); /* fixed levels of the CC work units */
    const int last = cc_unit_last(ctx, w, lev);
 
    atomic_store(&ctx->cc_target_workers, ctx->num_threads);
    const long long busy_ns0 = atomic_load(&ctx->cc_busy_ns);
    const long long nodes0 = atomic_load(&ctx->cc_nodes);
    atomic_store(&ctx->cc_next, cc_pack(w, lev, beg * ctx->cc_R[lev]));
    atomic_store(&ctx->cc_round_active, 1);
 
    double cc_start = get_time_sec();
    ctx->round_open = 1;
    ctx->round_w = w;
    ctx->round_beg = beg;
    ctx->round_end = end;
    ctx->round_R = ctx->cc_R[lev];
    ctx->round_threads = ctx->spawned;
    ctx->round_t0 = cc_start;
    ctx->round_est = t_cc_work;
 
    if (ctx->p->debug & DBG_TIMING) {
      const char *note = (certified || atomic_load(&ctx->cc_exact)) ? "" : " (lower weights not scanned)";
      char units[64] = "";
      if (lev > 0) snprintf(units, sizeof(units), ", %d work units per start column", ctx->cc_R[lev]);
      if (ctx->p->start_num > 0) {
        fprintf(stderr, "# searching w=%d with %d CC threads, start list (%d columns%s)%s%s\n",
                w, ctx->spawned, ctx->p->start_num, ctx->p->start_scan ? ", start_scan=1" : "", units, note);
      } else {
        fprintf(stderr, "# searching w=%d with %d CC threads, columns [%d, %d]%s%s\n",
                w, ctx->spawned, beg, end, units, note);
      }
    }
 
    bool round_completed = false;
    while (!atomic_load(&ctx->stop_flag)) {
      const double t = get_time_sec();
      if (ctx->timeout > 0.0 && (t - ctx->start_time >= ctx->timeout)) {
        ctx_stop(ctx, STOP_TIMEOUT);
        break;
      }
      if (cc_unpack_idx(atomic_load(&ctx->cc_next)) > last && atomic_load(&ctx->cc_active_workers) == 0) {
        round_completed = true;
        break;
      }
      if (m2_can_spawn(ctx) && t - cc_start > CC_SPAWN_TIME) {
        char why[64];
        snprintf(why, sizeof(why), "running for %.3fs", t - cc_start);
        m2_spawn_rest(ctx, w, why);
      }
      ctx_status(ctx);
      usleep(100);
    }
 
    atomic_store(&ctx->cc_round_active, 0);
    if (round_completed) ctx->round_open = 0;
 
    double cc_dur = get_time_sec() - cc_start;
    if (w < MAX_W) { /* CC work of the round in thread-seconds, as measured by the workers */
      ctx->cc_time_per_weight[w] = 1e-9 * (double)(atomic_load(&ctx->cc_busy_ns) - busy_ns0);
      if ((ctx->p->debug & DBG_TIMING) && round_completed) {
        char est_str[48] = "unknown";
        if (t_cc_work > 0.0) snprintf(est_str, sizeof(est_str), "%.3g thread-s", t_cc_work);
        fprintf(stderr, "# CC round w=%d work: %.3g thread-s measured, %s predicted (%.3fs on %d threads), %lld "
                "nodes\n", w, ctx->cc_time_per_weight[w], est_str, cc_dur, ctx->spawned,
                atomic_load(&ctx->cc_nodes) - nodes0);
      }
    }
 
    int cw_found = atomic_load(&ctx->cc_found_weight);
    if (cw_found > 0 && certified && cw_found == w && !atomic_load(&ctx->cc_exact)) {
      atomic_store(&ctx->cc_exact, true);
    }
    const bool exact = (cw_found > 0) && atomic_load(&ctx->cc_exact);
    if (exact) {
      atomic_store(&ctx->dmin, cw_found);
      atomic_store(&ctx->dmax, cw_found);
      /* the distance is known: no rounds beyond w = cw_found + extra_w */
      w_limit = minint(w_limit, cw_found + extra_w);
      const int next_w = (ctx->even_weight_only && !(w & 1)) ? w + 2 : w + 1;
      const bool more_rounds = round_completed && (next_w <= w_limit);
      if (ctx->p->debug & DBG_PROGRESS) {
        if (w == cw_found && more_rounds) {
          fprintf(stderr,
                  "# CC round w=%d finished in %.3fs (%d CC threads): found min-weight codewords "
                  "(dmin=%d, continuing up to w=%d for dW=%d, total %lld cws)\n",
                  w, cc_dur, ctx->spawned, cw_found, w_limit, ctx->p->dW, ctx_num_cws(ctx));
        } else if (w == cw_found) {
          fprintf(stderr, "# CC found min-weight codeword: d=%d (using %d CC threads, total %lld cws)\n",
                  cw_found, ctx->spawned, ctx_num_cws(ctx));
        } else if (round_completed) {
          fprintf(stderr,
                  "# CC round w=%d finished in %.3fs (%d CC threads): extra dW round completed "
                  "(dmin=%d, total %lld cws)\n",
                  w, cc_dur, ctx->spawned, cw_found, ctx_num_cws(ctx));
        }
      }
      if (!more_rounds) {
        /* (a stop reason recorded during the round, e.g., the timeout, is kept) */
        ctx_stop(ctx, collecting ? STOP_COLLECTED : STOP_CC_FOUND);
        break;
      }
    } else if (cw_found > 0) {
      /* noscan=1 without dmin=wmax: lower weights not scanned, the codeword only gives an upper bound */
      if (ctx->p->debug & DBG_PROGRESS) {
        fprintf(stderr,
                "# CC round w=%d finished in %.3fs (%d CC threads): found codeword of weight %d -> dmax=%d "
                "(dmin=%d not certified: lower weights not scanned)\n",
                w, cc_dur, ctx->spawned, cw_found, atomic_load(&ctx->dmax), atomic_load(&ctx->dmin));
      }
      ctx_stop(ctx, STOP_CC_FOUND);
      break;
    } else {
      if (!round_completed) {
        break;
      }
      if (certified) {
        /* Weight w analyzed without success */
        const int new_dmin = ctx->even_weight_only ? w + 2 : w + 1;
        atomic_store(&ctx->dmin, new_dmin);
        if (ctx->p->debug & DBG_PROGRESS) {
          fprintf(stderr, "# CC w=%d completed in %.3fs (%d CC threads): no codewords found -> dmin=%d\n",
                  w, cc_dur, ctx->spawned, new_dmin);
        }
      } else if (ctx->p->debug & DBG_PROGRESS) {
        fprintf(stderr, "# CC w=%d completed in %.3fs (%d CC threads): no codewords found "
                "(dmin=%d not raised: lower weights not scanned)\n",
                w, cc_dur, ctx->spawned, atomic_load(&ctx->dmin));
      }
    }
  }
}
 
/* Method 3: RW still runs (RW threads, RW steps not all completed, RW not stopped by min_hits) */
static inline bool m3_rw_running(distfork_ctx_t * const ctx) {
  return ctx->rw_threads > 0 && !atomic_load(&ctx->rw_stop_flag) &&
         atomic_load(&ctx->rw_steps_completed) < ctx->total_rw_steps;
}
 
/* Method 3: RW thread time per step, measured continuously by the RW workers (before the first step is completed, a
 * step takes at least the time elapsed since the start) */
static double m3_rw_step_time(distfork_ctx_t * const ctx) {
  const long n = atomic_load(&ctx->rw_timed_steps);
  if (n > 0) return 1e-9 * (double)atomic_load(&ctx->rw_busy_ns) / (double)n;
  const double elapsed = get_time_sec() - ctx->start_time;
  return (elapsed > 5e-5) ? elapsed : 5e-5;
}
 
/* Method 3: the largest CC weight needed: t-1 once the target t is known, i.e., the upper bound dmax or the stop
 * target dstop if smaller (with outC, t+dW: the rounds w = t..t+dW enumerate the codewords to export), otherwise wmax
 * (or n); at most MAX_W-2.  (`dexp` only affects the thread split, see m3_plan_cc_threads().) */
static int m3_cc_target_w(const distfork_ctx_t * const ctx, const int cur_dmin, const int cur_dmax) {
  const params_t * const p = ctx->p;
  const int t = ctx_target(ctx, cur_dmax);
  int target;
  if (t > 0) {
    if (p->outC && (p->dW > 0 || cur_dmin >= t)) {
      target = t + (p->dW > 0 ? p->dW : 0);
    } else {
      target = t - 1;
    }
  } else {
    target = (p->wmax > 0) ? p->wmax : p->spaH->cols;
  }
  const int max_allowed_w = (p->wmax > 0) ? minint(p->wmax, MAX_W - 2) : (MAX_W - 2);
  int res = minint(target, max_allowed_w);
  if (ctx->even_weight_only && (res & 1)) res--;
  return res;
}
 
#define M3_DEADLINE_FRACTION 0.75 /* with a timeout, a round gets enough CC threads to finish in 3/4 of the time left */
#define M3_OVERRUN_FACTOR 1.5     /* a round taking 1.5 times (+20 ms) longer than predicted gets all threads */
 
/* Method 3: number of CC threads for the round at weight w with the estimated work `est` (thread-seconds, 0 if not
 * known) and the growth factor `growth` (see cc_work_estimate()); the other threads run RW.  Returns 0 if CC is
 * paused (`dexp`, below). */
static int m3_plan_cc_threads(distfork_ctx_t * const ctx, const int w, const double est, const double growth,
                              const double remaining_time) {
  const int nthr = ctx->num_threads;
  if (nthr <= 1) return 1;
  /* RW finished, or all RW steps claimed: the RW workers join the round anyway */
  if (!m3_rw_running(ctx) || atomic_load(&ctx->rw_steps_started) >= ctx->total_rw_steps) return nthr;
  const int cur_dmin = atomic_load(&ctx->dmin);
  const int cur_dmax = atomic_load(&ctx->dmax);
  /* The round w = t-1 certifies dmin = t (t: dmax, or the stop target dstop if smaller; later rounds only collect
   * codewords).  RW cannot improve this result, it could only find a codeword of weight w before CC does: all threads
   * run CC. */
  const int t = ctx_target(ctx, cur_dmax);
  const int cert_w = ctx->even_weight_only ? (t - ((t & 1) ? 1 : 2)) : (t - 1);
  if (t > 0 && (w >= cert_w || cur_dmin >= t)) return nthr;
  /* `dexp` hint: as long as RW has not found any codeword, CC rounds w > dexp run only on the threads which cannot
   * run RW (none: CC pauses) */
  if (cur_dmax == 0 && ctx->dexp > 0 && w > ctx->dexp) return nthr - ctx->rw_threads;
 
  int n_cc;
  if (!(est > 0.0)) {
    n_cc = nthr / 2; /* the first round at a supplied dmin: its work is not known */
  } else if (est < 0.005) {
    n_cc = (nthr >= 4) ? 2 : 1;
  } else {
    /* Split by the remaining work in thread-seconds: the CC rounds w, w+1, w+2 up to the horizon (dmax-1, or dexp
     * before RW has found a codeword), and the remaining RW steps (at most 2000 once dmax is known), but not more RW
     * work than the RW threads can do before the timeout */
    int horizon = m3_cc_target_w(ctx, cur_dmin, cur_dmax);
    if (cur_dmax == 0 && ctx->dexp > 0 && ctx->dexp < horizon) horizon = ctx->dexp;
    double t_cc = est, t_k = est;
    const int step = ctx->even_weight_only ? 2 : 1;
    for (int k = w + step; k <= horizon && k <= w + 2 * step; k += step) {
      t_k *= growth;
      t_cc += t_k;
    }
    const long steps_rem = ctx->total_rw_steps - atomic_load(&ctx->rw_steps_completed);
    const long eff_steps = (cur_dmax > 0 && steps_rem > 2000) ? 2000 : steps_rem;
    double t_rw = (double)eff_steps * m3_rw_step_time(ctx);
    if (ctx->timeout > 0.0 && t_rw > remaining_time * ctx->rw_threads) t_rw = remaining_time * ctx->rw_threads;
    n_cc = (int)round((double)nthr * t_cc / (t_cc + t_rw));
    if (n_cc < 1) n_cc = 1;
    if (n_cc > nthr - 1) n_cc = nthr - 1; /* keep at least one RW thread */
  }
  /* timeout: enough CC threads to finish the round in time */
  if (ctx->timeout > 0.0 && est > 0.0 && remaining_time > 0.0) {
    const double n_min = ceil(est / (M3_DEADLINE_FRACTION * remaining_time));
    if (n_min > n_cc) n_cc = (n_min < nthr) ? (int)n_min : nthr;
  }
  /* only `rw_threads` workers can run RW (memory / steps limits): the other workers always run CC */
  if (nthr - n_cc > ctx->rw_threads) n_cc = nthr - ctx->rw_threads;
  return n_cc;
}
 
/* Method 3: wait while RW runs, until the upper bound dmax differs from `dmax0`, RW ends, or the run stops; `why`:
 * the reason why no CC round runs, for the status line */
static void m3_wait_rw(distfork_ctx_t * const ctx, const int dmax0, const char * const why) {
  ctx->cc_idle = why;
  while (!atomic_load(&ctx->stop_flag) && m3_rw_running(ctx) && atomic_load(&ctx->dmax) == dmax0) {
    if (ctx->timeout > 0.0 && (get_time_sec() - ctx->start_time >= ctx->timeout)) {
      ctx_stop(ctx, STOP_TIMEOUT);
      break;
    }
    ctx_status(ctx);
    usleep(1000);
  }
  ctx->cc_idle = NULL;
}
 
/* Method 3 coordinator */
static void run_method3_coordinator(distfork_ctx_t *ctx) {
  const int nvar = ctx->p->spaH->cols;
  const int nthr = ctx->num_threads;
  int w = ctx->p->noscan ? ctx->p->wmax : (ctx->p->dmin > 1 ? ctx->p->dmin : 1);
 
  int init_dmax = atomic_load(&ctx->dmax);
  int init_dmin = atomic_load(&ctx->dmin);
  if (init_dmax > 0 && init_dmin >= init_dmax && !ctx->p->outC) {
    ctx_stop(ctx, STOP_BOUNDS);
    return;
  }
 
  /* No RW probe: the RW step time is measured continuously by the RW workers (m3_rw_step_time()).  CC is never ended
   * while RW still runs: if no CC round can be started (CC done up to wmax, a round predicted to exceed the timeout,
   * or CC paused by `dexp`), the coordinator waits for RW and re-plans when dmax changes or RW ends. */
  int wait_msg_w = 0; /* weight w for which a waiting message was printed (once per round) */
  while (!atomic_load(&ctx->stop_flag)) {
    if (ctx->even_weight_only && (w & 1)) {
      if (w == atomic_load(&ctx->dmin)) atomic_store(&ctx->dmin, w + 1);
      w++;
      continue;
    }
    double now = get_time_sec();
    double remaining_time = (ctx->timeout > 0.0) ? (ctx->timeout - (now - ctx->start_time)) : 1e9;
    if (ctx->timeout > 0.0 && remaining_time <= 0.0) {
      ctx_stop(ctx, STOP_TIMEOUT);
      break;
    }
 
    int cur_dmax = atomic_load(&ctx->dmax);
    int cur_dmin = atomic_load(&ctx->dmin);
    const bool rw_on = m3_rw_running(ctx);
 
    /* Target cluster size for CC */
    const int target_cc_w = m3_cc_target_w(ctx, cur_dmin, cur_dmax);
 
    if (cur_dmax > 0 && cur_dmin >= cur_dmax && w > target_cc_w) {
      /* Bracketing converged and all requested dW rounds completed */
      atomic_store(&ctx->dmin, cur_dmax);
      ctx_stop(ctx, (ctx->p->outC && w > cur_dmax) ? STOP_COLLECTED : STOP_BOUNDS);
      break;
    }
    if (ctx->p->dstop > 0 && cur_dmin >= ctx->p->dstop && w > target_cc_w) {
      /* The stop target dstop reached (with outC, after the rounds w = dstop..dstop+dW which collect codewords) */
      ctx_stop(ctx, STOP_DSTOP);
      break;
    }
 
    if (w > target_cc_w) {
      /* CC done up to wmax: only RW can still lower dmax */
      if (!rw_on) {
        const int max_allowed_w = (ctx->p->wmax > 0) ? minint(ctx->p->wmax, MAX_W - 2) : (MAX_W - 2);
        ctx->stop_w = (ctx->even_weight_only && target_cc_w + 1 == max_allowed_w) ? max_allowed_w : target_cc_w;
        ctx_stop(ctx, STOP_CC_DONE);
        break;
      }
      if ((ctx->p->debug & DBG_TIMING) && wait_msg_w != w) {
        wait_msg_w = w;
        fprintf(stderr, "# CC done up to w=%d, waiting for RW (bounds [%d, %d])\n", target_cc_w, cur_dmin, cur_dmax);
      }
      m3_wait_rw(ctx, cur_dmax, "done up to wmax, waiting for RW");
      continue;
    }
 
    /* Estimated CC work of the round (thread-seconds).  A round predicted not to finish in time on all threads is not
     * started while RW runs; once RW has ended, the run ends.  Without RW (steps=0), as in method 2, the round is
     * started anyway: it can still end early with a codeword of weight w (exact if w = dmin). */
    double growth;
    const double t_cc_est = cc_work_estimate(ctx, w, &growth);
    if (ctx->timeout > 0.0 && t_cc_est / nthr > remaining_time) {
      if (rw_on) {
        if ((ctx->p->debug & DBG_TIMING) && wait_msg_w != w) {
          wait_msg_w = w;
          fprintf(stderr, "# CC for w=%d (est %.2fs) exceeds remaining timeout %.2fs, devoting %d threads to RW\n",
                  w, t_cc_est / nthr, remaining_time, ctx->rw_threads);
        }
        m3_wait_rw(ctx, cur_dmax, "round predicted to exceed the timeout, waiting for RW");
        continue;
      }
      if (ctx->rw_threads > 0) {
        if (ctx->p->debug & DBG_PROGRESS) {
          fprintf(stderr, "# CC for w=%d (est %.2fs) exceeds remaining timeout %.2fs, RW ended: terminating early "
                  "(dmin=%d)\n", w, t_cc_est / nthr, remaining_time, cur_dmin);
        }
        ctx->stop_w = w;
        ctx_stop(ctx, STOP_CC_SLOW);
        break;
      }
      if ((ctx->p->debug & DBG_PROGRESS) && wait_msg_w != w) {
        wait_msg_w = w;
        fprintf(stderr, "# CC for w=%d (est %.2fs) may not finish in the remaining time %.2fs, searching anyway "
                "(dmin=%d)\n", w, t_cc_est / nthr, remaining_time, cur_dmin);
      }
    }
 
    /* Calculate thread balancing */
    int n_cc = m3_plan_cc_threads(ctx, w, t_cc_est, growth, remaining_time);
    if (n_cc <= 0) {
      if ((ctx->p->debug & DBG_TIMING) && wait_msg_w != w) {
        wait_msg_w = w;
        fprintf(stderr, "# CC paused before w=%d > dexp=%d until RW finds a codeword (or ends)\n", w, ctx->dexp);
      }
      m3_wait_rw(ctx, cur_dmax, "paused (w > dexp) until RW finds a codeword");
      continue;
    }
    const int n_cc0 = n_cc;
    int n_rw = nthr - n_cc;
    const long steps_rem = rw_on ? ctx->total_rw_steps - atomic_load(&ctx->rw_steps_completed) : 0;
 
    const int beg = cc_round_beg(ctx->p);
    const int end = cc_round_end(ctx->p, w);
    /* all weights below w analyzed (or a supplied dmin): a codeword of weight w found is exact */
    const bool certified = (w == atomic_load(&ctx->dmin));
    const int lev = cc_round_levels(ctx, w, beg, end, t_cc_est); /* fixed levels of the CC work units */
    const int last = cc_unit_last(ctx, w, lev);
 
    atomic_store(&ctx->cc_target_workers, n_cc);
    const long long busy_ns0 = atomic_load(&ctx->cc_busy_ns);
    const long long nodes0 = atomic_load(&ctx->cc_nodes);
    atomic_store(&ctx->cc_next, cc_pack(w, lev, beg * ctx->cc_R[lev]));
    atomic_store(&ctx->cc_round_active, 1);
 
    if (ctx->p->debug & DBG_TIMING) {
      char est_str[64] = "unknown";
      char rw_str[64] = "";
      char units[64] = "";
      if (t_cc_est > 0.0) snprintf(est_str, sizeof(est_str), "%.3g thread-s", t_cc_est);
      if (rw_on) snprintf(rw_str, sizeof(rw_str), ", RW step %.3g s", m3_rw_step_time(ctx));
      if (lev > 0) snprintf(units, sizeof(units), ", %d work units per start column", ctx->cc_R[lev]);
      fprintf(stderr,
              "# CC round w=%d started: %d CC threads, %d RW threads "
              "(bounds [%d, %d], rem_rw=%ld, rem_time=%.2fs, est %s%s%s)\n",
              w, n_cc, n_rw, cur_dmin, cur_dmax, steps_rem, remaining_time, est_str, rw_str, units);
    }
 
    double cc_start = get_time_sec();
    bool round_completed = false;
    int dmax_seen = cur_dmax;
    unsigned int poll = 0;
    ctx->round_open = 1;
    ctx->round_w = w;
    ctx->round_beg = beg;
    ctx->round_end = end;
    ctx->round_R = ctx->cc_R[lev];
    ctx->round_threads = n_cc;
    ctx->round_t0 = cc_start;
    ctx->round_est = t_cc_est;
 
    while (!atomic_load(&ctx->stop_flag)) {
      const double t = get_time_sec();
      if (ctx->timeout > 0.0 && (t - ctx->start_time >= ctx->timeout)) {
        ctx_stop(ctx, STOP_TIMEOUT);
        break;
      }
      if (cc_unpack_idx(atomic_load(&ctx->cc_next)) > last && atomic_load(&ctx->cc_active_workers) == 0) {
        round_completed = true;
        break;
      }
      /* Re-plan during the round (only adding CC threads, workers do not leave a round): when RW finds a new upper
       * bound (e.g., this round now certifies dmin = dmax), and when the round takes much longer than predicted */
      if ((++poll & 15) == 0 && n_cc < nthr) {
        int n_new = n_cc;
        const char *why = "";
        const int d_now = atomic_load(&ctx->dmax);
        if (d_now != dmax_seen) {
          dmax_seen = d_now;
          const double rem = (ctx->timeout > 0.0) ? ctx->timeout - (t - ctx->start_time) : 1e9;
          n_new = m3_plan_cc_threads(ctx, w, t_cc_est, growth, rem);
          why = "new upper bound";
        }
        const bool paused = (d_now == 0 && ctx->dexp > 0 && w > ctx->dexp && m3_rw_running(ctx));
        if (!paused && t_cc_est > 0.0 && t - cc_start > M3_OVERRUN_FACTOR * t_cc_est / n_cc0 + 0.02) {
          n_new = nthr;
          why = "longer than predicted";
        }
        if (n_new > n_cc) {
          if (ctx->p->debug & DBG_TIMING) {
            fprintf(stderr, "# CC round w=%d: %d -> %d CC threads (%s, bounds [%d, %d], %.3fs into the round)\n",
                    w, n_cc, n_new, why, atomic_load(&ctx->dmin), d_now, t - cc_start);
          }
          n_cc = n_new;
          n_rw = nthr - n_cc;
          ctx->round_threads = n_cc;
          atomic_store(&ctx->cc_target_workers, n_cc);
        }
      }
      ctx_status(ctx);
      usleep(100);
    }
 
    atomic_store(&ctx->cc_round_active, 0);
    if (round_completed) ctx->round_open = 0;
 
    double cc_dur = get_time_sec() - cc_start;
    if (w < MAX_W) {
      /* CC work of the round in thread-seconds, as measured by the workers: the number of CC threads may exceed
       * n_cc (CC-only workers, and RW workers that ran out of RW steps join the round) */
      ctx->cc_time_per_weight[w] = 1e-9 * (double)(atomic_load(&ctx->cc_busy_ns) - busy_ns0);
      if ((ctx->p->debug & DBG_TIMING) && round_completed) {
        char est_str[48] = "unknown";
        if (t_cc_est > 0.0) snprintf(est_str, sizeof(est_str), "%.3g thread-s", t_cc_est);
        fprintf(stderr, "# CC round w=%d work: %.3g thread-s measured, %s predicted (%.3fs, %d CC threads at the "
                "end), %lld nodes\n", w, ctx->cc_time_per_weight[w], est_str, cc_dur, n_cc,
                atomic_load(&ctx->cc_nodes) - nodes0);
      }
    }
 
    int cw_found = atomic_load(&ctx->cc_found_weight);
    if (cw_found > 0 && certified && cw_found == w && !atomic_load(&ctx->cc_exact)) {
      atomic_store(&ctx->cc_exact, true);
    }
    const bool exact = (cw_found > 0) && atomic_load(&ctx->cc_exact);
    if (exact) {
      atomic_store(&ctx->dmin, cw_found);
      atomic_store(&ctx->dmax, cw_found);
 
      int max_w_lim = minint(ctx->p->wmax > 0 ? ctx->p->wmax : nvar, cw_found + ctx->p->dW);
      if (ctx->p->dstop > 0 && ctx->p->dstop < cw_found) {
        /* a stop target dstop below the distance: the codewords of weight up to dstop + dW only (as in method 2) */
        max_w_lim = minint(max_w_lim, ctx->p->dstop + (ctx->p->dW > 0 ? ctx->p->dW : 0));
      }
      if (ctx->even_weight_only && (max_w_lim & 1)) max_w_lim--;
      if (ctx->p->outC && ctx->p->dW > 0 && w < max_w_lim) {
        if (ctx->p->debug & DBG_PROGRESS) {
          if (w == cw_found) {
            fprintf(stderr,
                    "# CC round w=%d finished in %.3fs (%d CC threads, %d RW threads): found codewords "
                    "(dmin=%d, continuing up to w=%d for dW=%d, total %lld cws)\n",
                    w, cc_dur, n_cc, n_rw, cw_found, max_w_lim, ctx->p->dW, ctx_num_cws(ctx));
          } else if (round_completed) {
            fprintf(stderr,
                    "# CC round w=%d finished in %.3fs (%d CC threads, %d RW threads): "
                    "extra dW round completed (dmin=%d, total %lld cws)\n",
                    w, cc_dur, n_cc, n_rw, cw_found, ctx_num_cws(ctx));
          }
        }
      } else {
        if (ctx->p->debug & DBG_PROGRESS) {
          if (w > cw_found) {
            if (round_completed) {
              fprintf(stderr,
                      "# CC round w=%d finished in %.3fs (%d CC threads, %d RW threads): "
                      "extra dW round completed (dmin=%d, total %lld cws)\n",
                      w, cc_dur, n_cc, n_rw, cw_found, ctx_num_cws(ctx));
            }
          } else {
            fprintf(stderr, "# CC found min-weight codeword: d=%d (using %d CC threads, total %lld cws)\n",
                    cw_found, n_cc, ctx_num_cws(ctx));
          }
        }
        /* (a stop reason recorded during the round, e.g., the timeout, is kept) */
        ctx_stop(ctx, (ctx->p->outC || ctx->p->maxC > 0) ? STOP_COLLECTED : STOP_CC_FOUND);
        break;
      }
    } else if (cur_dmax > 0 && cur_dmin >= cur_dmax) {
      /* Extra dW round completed */
      if (round_completed && (ctx->p->debug & DBG_PROGRESS)) {
        fprintf(stderr,
                "# CC round w=%d finished in %.3fs (%d CC threads, %d RW threads): "
                "extra dW round completed (dmin=%d, total %lld cws)\n",
                w, cc_dur, n_cc, n_rw, cur_dmin, ctx_num_cws(ctx));
      }
      if (!round_completed) {
        break;
      }
    } else {
      if (!round_completed) {
        break;
      }
 
      /* Weight w analyzed without success */
      int new_dmin = ctx->even_weight_only ? w + 2 : w + 1;
      atomic_store(&ctx->dmin, new_dmin);
      if (ctx->p->debug & DBG_PROGRESS) {
        fprintf(stderr, "# CC round w=%d finished in %.3fs (%d CC threads, %d RW threads): no codewords -> dmin=%d\n",
                w, cc_dur, n_cc, n_rw, new_dmin);
      }
 
      cur_dmax = atomic_load(&ctx->dmax);
      if (cur_dmax > 0 && new_dmin >= cur_dmax) {
        atomic_store(&ctx->dmin, cur_dmax);
        /* with outC, the CC rounds w = dmax..dmax+dW enumerate all codewords to export (RW may have missed some) */
        const int w_last = ctx->p->outC ? m3_cc_target_w(ctx, cur_dmax, cur_dmax) : 0;
        if (w_last > w) {
          if (ctx->p->debug & DBG_PROGRESS) {
            fprintf(stderr, "# bracketing bounds coincide: dmin = dmax = %d (continuing up to w=%d for dW=%d to export "
                    "codewords)\n", cur_dmax, w_last, ctx->p->dW > 0 ? ctx->p->dW : 0);
          }
        } else {
          ctx_stop(ctx, STOP_BOUNDS);
          if (ctx->p->debug & DBG_PROGRESS) {
            fprintf(stderr, "# bracketing bounds coincide: dmin = dmax = %d\n", cur_dmax);
          }
          break;
        }
      } else if (ctx->p->dstop > 0 && new_dmin >= ctx->p->dstop) {
        /* the stop target dstop reached: no codeword of weight < dstop (RW may still find heavier codewords, but they
         * are not needed); with outC, the CC rounds w = dstop..dstop+dW collect codewords first */
        const int w_last = ctx->p->outC ? m3_cc_target_w(ctx, new_dmin, cur_dmax) : 0;
        if (w_last > w) {
          if (ctx->p->debug & DBG_PROGRESS) {
            fprintf(stderr, "# lower bound dmin=%d reached dstop=%d (continuing up to w=%d for dW=%d to export "
                    "codewords)\n", new_dmin, ctx->p->dstop, w_last, ctx->p->dW > 0 ? ctx->p->dW : 0);
          }
        } else {
          ctx_stop(ctx, STOP_DSTOP);
          if (ctx->p->debug & DBG_PROGRESS) {
            fprintf(stderr, "# lower bound dmin=%d reached dstop=%d\n", new_dmin, ctx->p->dstop);
          }
          break;
        }
      }
    }
 
    w++;
  }
}
 
/* Default number of worker threads: the logical CPUs the process may run on, i.e., its CPU affinity mask (`taskset`,
 * cpusets) on Linux, otherwise the online CPUs (4 if not known).  CPU quotas (e.g., `docker --cpus`) are not seen. */
static int available_cpus(void) {
#ifdef __linux__
  cpu_set_t set;
  CPU_ZERO(&set);
  if (sched_getaffinity(0, sizeof(set), &set) == 0) { /* fails with more than CPU_SETSIZE (1024) CPUs */
    const int cnt = CPU_COUNT(&set);
    if (cnt > 0) return cnt;
  }
#endif
  const long nprocs = sysconf(_SC_NPROCESSORS_ONLN);
  return (nprocs > 0) ? (int)nprocs : 4;
}
 
/* The physical memory in bytes (0 if not known) */
static unsigned long long physical_memory_bytes(void) {
#if defined(_SC_PHYS_PAGES) && defined(_SC_PAGESIZE)
  const long pages = sysconf(_SC_PHYS_PAGES), psize = sysconf(_SC_PAGESIZE);
  if (pages > 0 && psize > 0) return (unsigned long long)pages * (unsigned long long)psize;
#endif
  return 0;
}
 
/* Maximum number of RW workers for thread throttling (all workers in method 1, the RW share in method 3):
 * - dense memory: a full-matrix RW worker holds dense copies of H (m x n) and of its transpose (n x m), a worker
 *   with ksub > 0 only the sampled subspace matrix (ksub x n).  Above 15 MB per worker, RW is limited by the memory
 *   bandwidth: the throughput peaks at 32-64 RW workers whatever the matrix size (e.g., 32 with 49 MB per worker,
 *   32-64 with 0.43 GB).  Then at most 32 workers run RW, with their dense matrices in a quarter of the physical memory
 *   (but at least 2 workers; 1.5 GB if the physical memory is not known);
 * - small step counts: at most ceil(steps/10) RW workers. */
static int rw_thread_cap(const params_t * const p) {
  int cap = INT_MAX;
  const int nvar = p->nvar;
  const int nrows = p->spaH ? p->spaH->rows : 0;
  if (nrows > 0 && nvar > 0) {
    const size_t words_n = ((size_t)nvar + 63) / 64;
    size_t dense_bytes_per_thread;
    if (p->ksub > 0) {
      dense_bytes_per_thread = (size_t)minint(p->ksub, nrows) * words_n * 8;
    } else {
      dense_bytes_per_thread = (size_t)nrows * words_n * 8 + (size_t)nvar * (((size_t)nrows + 63) / 64) * 8;
    }
    if (dense_bytes_per_thread > 15ULL * 1024 * 1024) {
      const unsigned long long ram = physical_memory_bytes();
      const unsigned long long budget = ram ? ram / 4 : 1536ULL * 1024 * 1024;
      const unsigned long long max_mem_threads = budget / dense_bytes_per_thread;
      cap = (max_mem_threads < 2) ? 2 : ((max_mem_threads > 32) ? 32 : (int)max_mem_threads);
    }
  }
  if (p->steps > 0) {
    const int max_step_threads = (int)(((long)p->steps + 9) / 10);
    cap = minint(cap, max_step_threads);
  }
  return cap;
}
 
/* DBG_PROGRESS: the run plan (method, threads, RW steps, CC weights, timeout, rng seed) */
static void print_run_plan(const distfork_ctx_t * const ctx, const int init_dmax) {
  const params_t * const p = ctx->p;
  fprintf(stderr, "# run plan: method=%d (%s), %d threads", p->method,
          (p->method == 1) ? "RW" : ((p->method == 2) ? "CC" : "bracketing"), ctx->num_threads);
  if (ctx->spawn_initial < ctx->num_threads) fprintf(stderr, " (%d at first)", ctx->spawn_initial);
  if (p->method == 3) fprintf(stderr, " (RW on up to %d)", ctx->rw_threads);
  if ((p->method & 1) && ctx->rw_threads > 0)
    fprintf(stderr, ", steps=%ld, min_hits=%d, ksub=%d", ctx->total_rw_steps, p->min_hits, p->ksub);
  if (p->method & 2) {
    fprintf(stderr, ", CC from w=%d", p->noscan ? p->wmax : (p->dmin > 1 ? p->dmin : 1));
    if (p->wmax > 0) fprintf(stderr, " to wmax=%d", p->wmax);
    if (p->smax == 0) fprintf(stderr, " (cc_prune=%d)", p->cc_prune);
    if (p->method == 3 && ctx->dexp > 0) fprintf(stderr, ", dexp=%d", ctx->dexp);
  }
  if (init_dmax > 0) fprintf(stderr, ", dmax=%d", init_dmax);
  if (p->dstop > 0) fprintf(stderr, ", dstop=%d", p->dstop);
  if (ctx->timeout > 0.0) fprintf(stderr, ", timeout=%gs", ctx->timeout);
  else fprintf(stderr, ", no timeout");
  fprintf(stderr, ", seed=%d\n", p->seed);
}
 
/* DBG_SUMMARY: why the run ended (see ctx_stop()) */
static void print_stop_reason(distfork_ctx_t * const ctx, const int final_dmax, const bool cc_exact) {
  params_t * const p = ctx->p;
  const int cc_found = atomic_load(&ctx->cc_found_weight);
  const bool rw_ran = (p->method & 1) && ctx->rw_threads > 0;
  fprintf(stderr, "# stopped after %.3fs: ", get_time_sec() - ctx->start_time);
  switch (atomic_load(&ctx->stop_reason)) {
  case STOP_BOUNDS:
    fprintf(stderr, "bounds coincide: dmin = dmax = %d", final_dmax);
    break;
  case STOP_CC_FOUND:
    if (cc_exact)
      fprintf(stderr, "CC found min-weight codeword: d=%d", cc_found);
    else
      fprintf(stderr, "CC found a codeword of weight %d (upper bound only: lower weights not scanned)", cc_found);
    break;
  case STOP_RW_CONVERGED: {
    cw_hit_stats_t st;
    codeword_hit_stats(p, &st);
    fprintf(stderr, "RW convergence reached (<n>=%.2f >= min_hits=%d after %ld RW steps)", st.avg, p->min_hits,
            atomic_load(&ctx->rw_steps_completed));
    break;
  }
  case STOP_RW_STEPS:
    fprintf(stderr, "all %ld RW steps done", ctx->total_rw_steps);
    break;
  case STOP_TIMEOUT:
    fprintf(stderr, "timeout=%gs reached", ctx->timeout);
    if (ctx->round_open)
      fprintf(stderr, " during CC round w=%d (%.1f%% of start columns claimed)", ctx->round_w, cc_round_claimed(ctx));
    if (rw_ran)
      fprintf(stderr, ", %ld of %ld RW steps done", atomic_load(&ctx->rw_steps_completed), ctx->total_rw_steps);
    break;
  case STOP_CC_DONE:
    fprintf(stderr, "CC done up to w=%d", ctx->stop_w);
    if (p->method == 3 && rw_ran) fprintf(stderr, ", RW ended (%s)", rw_end_state(ctx));
    break;
  case STOP_CC_SLOW:
    fprintf(stderr, "CC round w=%d predicted not to finish before the timeout, RW ended (%s)", ctx->stop_w,
            rw_end_state(ctx));
    break;
  case STOP_WMIN:
    fprintf(stderr, "found a codeword of weight %d <= wmin=%d", final_dmax, p->wmin);
    break;
  case STOP_MAXC:
    fprintf(stderr, "maxC=%lld codewords collected", p->maxC);
    break;
  case STOP_COLLECTED: {
    /* the collection rounds w = d..d+dW (with outC), but at most up to dstop+dW (a stop target dstop < d) and wmax */
    int w_hi = final_dmax;
    if (p->outC && p->dW > 0) {
      w_hi = final_dmax + p->dW;
      if (p->dstop > 0 && p->dstop < final_dmax) w_hi = minint(w_hi, p->dstop + p->dW);
      if (p->wmax > 0) w_hi = minint(w_hi, p->wmax);
    }
    fprintf(stderr, "CC enumerated all codewords of weight %d", final_dmax);
    if (w_hi > final_dmax) fprintf(stderr, "..%d", w_hi);
    if (p->outC) fprintf(stderr, " (for outC)");
    break;
  }
  case STOP_DSTOP:
    fprintf(stderr, "lower bound dmin=%d reached dstop=%d", atomic_load(&ctx->dmin), p->dstop);
    break;
  default:
    fprintf(stderr, "search ended");
  }
  fprintf(stderr, "\n");
}
 
typedef struct {
  uint64_t h;
  int c;
} cc_hash_col_t;
 
static int cc_hash_col_cmp(const void *a, const void *b) {
  const cc_hash_col_t *x = a, *y = b;
  if (x->h != y->h) return (x->h < y->h) ? -1 : 1;
  return (x->c > y->c) - (x->c < y->c);
}
 
/* The CC column-syndrome hash of H (mHT: its transpose with m rows), see cc_hash_t; NULL if H has no nonzero column */
static cc_hash_t *cc_hash_new(const csr_t * const mHT, const int m) {
  const int n = mHT->rows;
  cc_hash_t * const hs = calloc(1, sizeof(cc_hash_t));
  cc_hash_col_t * const tmp = malloc((size_t)(n > 0 ? n : 1) * sizeof(cc_hash_col_t));
  if (!hs || !tmp) ERROR("memory allocation");
  hs->zr = malloc((size_t)(m > 0 ? m : 1) * sizeof(uint64_t));
  hs->hc = calloc((size_t)(n > 0 ? n : 1), sizeof(uint64_t));
  if (!hs->zr || !hs->hc) ERROR("memory allocation");
  uint64_t state = 0x243f6a8885a308d3ULL; /* fixed: the results do not depend on the hash */
  for (int r = 0; r < m; r++) hs->zr[r] = splitmix64(&state);
  int nz = 0;
  for (int c = 0; c < n; c++) {
    if (mHT->p[c + 1] == mHT->p[c]) continue;
    uint64_t h = 0;
    for (int k = mHT->p[c]; k < mHT->p[c + 1]; k++) h ^= hs->zr[mHT->i[k]];
    hs->hc[c] = h;
    tmp[nz].h = h;
    tmp[nz++].c = c;
  }
  if (nz == 0) {
    free(tmp);
    free(hs->zr);
    free(hs->hc);
    free(hs);
    return NULL;
  }
  qsort(tmp, nz, sizeof(cc_hash_col_t), cc_hash_col_cmp);
  int nkeys = 0;
  for (int k = 0; k < nz; k++) nkeys += (k == 0 || tmp[k].h != tmp[k - 1].h);
  int lg = 0; /* ceil(log2(nkeys)) */
  while ((1LL << lg) < nkeys) lg++;
  const int tbits = lg + 1;             /* table: at most half full */
  const int fbits = (lg + 4 > 6) ? lg + 4 : 6; /* filter: 16 bits per key */
  hs->tshift = 64 - tbits;
  hs->tmask = (int)((1ULL << tbits) - 1);
  hs->tab = calloc((size_t)1 << tbits, sizeof(cc_hash_slot_t));
  hs->bmask = ((uint64_t)1 << fbits) - 1;
  hs->bits = calloc((size_t)1 << (fbits - 6), sizeof(uint64_t));
  hs->cols = malloc((size_t)nz * sizeof(int));
  if (!hs->tab || !hs->bits || !hs->cols) ERROR("memory allocation");
  for (int k = 0; k < nz; k++) hs->cols[k] = tmp[k].c;
  for (int k = 0; k < nz;) {
    const uint64_t key = tmp[k].h;
    int k2 = k + 1;
    while (k2 < nz && tmp[k2].h == key) k2++;
    int j = (int)(key >> hs->tshift);
    while (hs->tab[j].cnt > 0) j = (j + 1) & hs->tmask;
    hs->tab[j].key = key;
    hs->tab[j].beg = k;
    hs->tab[j].cnt = k2 - k;
    const uint64_t fb = key & hs->bmask;
    hs->bits[fb >> 6] |= (uint64_t)1 << (fb & 63);
    k = k2;
  }
  hs->ncols = nz;
  hs->nkeys = nkeys;
  free(tmp);
  return hs;
}
 
static void cc_hash_free(cc_hash_t * const hs) {
  if (!hs) return;
  free(hs->zr);
  free(hs->hc);
  free(hs->bits);
  free(hs->tab);
  free(hs->cols);
  free(hs);
}
 
int main(int argc, char **argv) {
  params_t * const p = &prm;
 
  var_init(argc, argv, p);
 
  if (p->finC) {
    nzlist_read(p->finC, p);
  }
 
  /* Check ksub feasibility before thread throttling (require m >= nu >= n - m) */
  if ((p->method & 1) && p->ksub > 0 && p->spaH) {
    int nrows = p->spaH->rows;
    int nvar = p->spaH->cols;
    if (2 * nrows < nvar) {
      if (p->debug & DBG_SUMMARY) {
        fprintf(stderr,
                "# Warning: ksub=%d requested, but m=%d < n-m=%d (<= nu); "
                "falling back to full-matrix RW (ksub=0)\n",
                p->ksub, nrows, nvar - nrows);
      }
      p->ksub = 0;
    }
  }
 
  /* Determine number of threads (default: all logical CPUs available to the process) */
  int requested_threads = p->threads;
  int num_threads = p->threads;
  if (num_threads <= 0) num_threads = available_cpus();
  /* Workers which may run RW (each allocates dense RW matrices): all workers in method 1, none in method 2.
   * In method 3, the RW share may be limited further while CC rounds can use all workers. */
  int rw_threads = 0;
 
  if (!p->nothrottle) {
    /* Thread throttling for small codes and large-memory matrices */
    int nvar = p->nvar;
    int nrows = p->spaH ? p->spaH->rows : 0;
    unsigned long long n_elements = (unsigned long long)nrows * (unsigned long long)nvar;
 
    /* 1. Small-code throttling by the RW work (methods 1 and 3): avoid the thread start-up overhead on tiny
     *    workloads.  The CC work is not known in advance: method 2 starts a few workers only, and the others once a
     *    CC round takes longer (see CC_SPAWN_INITIAL) */
    if (p->method & 1) {
      const double rw_work = (double)n_elements * (double)(p->steps > 0 ? p->steps : 1);
      if (nvar < 60 || (n_elements < 20000ULL && rw_work < 5.0e7)) {
        if (num_threads > 4) num_threads = 4;
      } else if ((nvar < 150 || n_elements < 100000ULL) && rw_work < 2.0e8) {
        if (num_threads > 16) num_threads = 16;
      }
    }
 
    /* 2. Large-matrix memory throttling and 3. steps throttling: these limit only the RW workers, i.e., all
     *    workers in method 1 and the RW share in method 3; CC (method 2, CC rounds in method 3) is not limited */
    if (p->method & 1) {
      const int rw_cap = rw_thread_cap(p);
      if (p->method == 1 && num_threads > rw_cap) num_threads = rw_cap;
      rw_threads = minint(num_threads, rw_cap);
    }
 
    if (requested_threads > 0 && num_threads != requested_threads && (p->debug & DBG_TIMING)) {
      fprintf(stderr, "# Note: throttled threads from %d to %d (n=%d, r=%d)\n",
              requested_threads, num_threads, nvar, nrows);
    }
  } else {
    if (p->method & 1) rw_threads = num_threads;
    if (p->debug & DBG_TIMING) {
      fprintf(stderr, "# Thread throttling disabled (nothrottle=1); running with %d threads\n", num_threads);
    }
  }
  if (p->method == 3 && p->steps == 0) {
    rw_threads = 0; /* no RW steps: pure CC through the bracketing coordinator, no RW matrices needed */
  }
  if (p->method == 3 && rw_threads < num_threads && (p->debug & DBG_TIMING)) {
    if (rw_threads == 0) {
      fprintf(stderr, "# Note: steps=0, no RW; CC rounds use all %d threads\n", num_threads);
    } else {
      fprintf(stderr, "# Note: RW limited to %d of %d threads (memory, steps); CC rounds can use all threads\n",
              rw_threads, num_threads);
    }
  }
 
  double timeout = (p->timeout > 0.0) ? p->timeout : 0.0;
 
  distfork_ctx_t ctx;
  memset(&ctx, 0, sizeof(ctx));
  ctx.p = p;
  ctx.num_threads = num_threads;
  ctx.rw_threads = rw_threads;
  ctx.timeout = timeout;
  ctx.start_time = get_time_sec();
  ctx.dexp = p->dexp;
  ctx.total_rw_steps = (p->steps >= 0) ? p->steps : 1;
 
  /* Initialize dmin and dmax */
  atomic_init(&ctx.dmin, p->dmin > 1 ? p->dmin : 1);
  int init_dmax = 0;
  if (p->dmax > 0) {
    init_dmax = p->dmax;
  }
  if (p->min_w != INT_MAX) {
    if (init_dmax == 0 || p->min_w < init_dmax) {
      init_dmax = p->min_w;
    }
  }
  atomic_init(&ctx.dmax, init_dmax);
 
  if (init_dmax > 0 && p->wmin > 0 && init_dmax <= p->wmin && !p->outC && p->maxC == 0) {
    if (p->debug & DBG_SUMMARY) {
      fprintf(stderr, "# stopped: known upper bound dmax=%d <= wmin=%d (no search)\n", init_dmax, p->wmin);
    }
    printf("%d %d 0\n", p->dmin > 1 ? p->dmin : 1, init_dmax);
    var_kill(p);
    return 0;
  }
 
  if (p->method == 3 && init_dmax > 0 && p->dmin > 1 && p->dmin >= init_dmax && !p->outC) {
    if (p->debug & DBG_SUMMARY) {
      fprintf(stderr, "# stopped: bounds coincide: dmin = dmax = %d (supplied, no search)\n", init_dmax);
    }
    printf("%d %d 0\n", p->dmin, init_dmax);
    var_kill(p);
    return 0;
  }
 
  /* the supplied dmin already reaches the stop target dstop (any method; with outC or maxC, the CC rounds which collect
   * codewords still run) */
  if (p->dstop > 0 && p->dmin >= p->dstop && !p->outC && p->maxC == 0) {
    if (p->debug & DBG_SUMMARY) {
      fprintf(stderr, "# stopped: lower bound dmin=%d >= dstop=%d (supplied, no search)\n", p->dmin, p->dstop);
    }
    printf("%d %d 0\n", (init_dmax > 0 && p->dmin >= init_dmax) ? init_dmax : p->dmin, init_dmax);
    var_kill(p);
    return 0;
  }
 
  atomic_init(&ctx.cc_found_weight, 0);
  atomic_init(&ctx.cc_exact, false);
  atomic_init(&ctx.stop_flag, false);
  atomic_init(&ctx.rw_stop_flag, false);
  atomic_init(&ctx.stop_reason, STOP_NONE);
  ctx.status_next = 1.0; /* DBG_STATUS: after 1, 2, 4, ... s, then every 60 s */
  ctx.status_dt = 1.0;
  atomic_init(&ctx.cw_feed_limit, codeword_feed_limit(p));
  atomic_init(&ctx.rw_steps_started, 0);
  atomic_init(&ctx.rw_steps_completed, 0);
  atomic_init(&ctx.rw_busy_ns, 0);
  atomic_init(&ctx.rw_timed_steps, 0);
  atomic_init(&ctx.rw_rank, 0);
  atomic_init(&ctx.rw_unif_steps, 0);
  atomic_init(&ctx.rw_pairs_l, 0);
  atomic_init(&ctx.rw_pairs_checked, 0);
  atomic_init(&ctx.rw_hits1, 0);
  atomic_init(&ctx.rw_hits2, 0);
  uint64_t * const rw_lhash = (p->pairs != 0 && (p->method & 1) && p->ksub <= 0) ? rw_lhash_init(p) : NULL;
  ctx.rw_lhash = rw_lhash;
  atomic_init(&ctx.cc_next, cc_pack(0, 0, 0));
  atomic_init(&ctx.cc_active_workers, 0);
  atomic_init(&ctx.cc_target_workers, 0);
  atomic_init(&ctx.cc_round_active, 0);
  atomic_init(&ctx.cc_busy_ns, 0);
  atomic_init(&ctx.cc_nodes, 0);
  atomic_init(&ctx.cw_red_cc, 0);
  atomic_init(&ctx.cw_red_rw, 0);
  atomic_init(&ctx.next_refresh_step, p->refresh > 0 ? p->refresh : 0);
 
  pthread_mutex_init(&ctx.cw_mutex, NULL);
  { /* ksub basis lock: writer-preferring (glibc), so that a waiting refresh (writer) is not starved by overlapping RW
     * batches (readers, which never hold the lock recursively) */
    pthread_rwlockattr_t rwattr;
    pthread_rwlockattr_init(&rwattr);
#ifdef __GLIBC__
    pthread_rwlockattr_setkind_np(&rwattr, PTHREAD_RWLOCK_PREFER_WRITER_NONRECURSIVE_NP);
#endif
    pthread_rwlock_init(&ctx.basis_rwlock, &rwattr);
    pthread_rwlockattr_destroy(&rwattr);
  }
 
  ctx.mHT_cc = csr_transpose(NULL, p->spaH);
  ctx.max_col_W = csr_max_row_wght(ctx.mHT_cc);
  /* the CC walk with bitset syndromes and the column-syndrome hash for its last two levels (see cc_walk(); not with
   * smax > 0, which needs the syndrome weights of the leaf level) */
  cc_hash_t * const cc_hash = ((p->method & 2) && p->smax == 0) ? cc_hash_new(ctx.mHT_cc, p->spaH->rows) : NULL;
  ctx.cc_hash = cc_hash;
  ctx.even_weight_only = ((p->method & 2) && p->smax == 0 && !p->noscan && p->spaH &&
                          csr_all_ones_in_rowspan(p->spaH));
  if (cc_hash && (p->debug & DBG_TIMING)) {
    fprintf(stderr, "# CC: %s, bitset syndromes, hash of %d nonzero columns (%d distinct) for the last two levels"
            " (%d table slots)%s\n", (p->cc_prune == 0) ? "walk of the first unsatisfied row (cc_prune=0)"
            : ((p->cc_prune == 1) ? "branch and bound (cc_prune=1)" : "branch and bound with the cover bounds "
               "(cc_prune=2)"), cc_hash->ncols, cc_hash->nkeys, cc_hash->tmask + 1,
            ctx.even_weight_only ? ", even weights only (1 in rowspan(H))" : "");
  }
  /* CC work units (see `cc_pack()`): radix B = the maximum row weight of H, B^L units per start index */
  ctx.cc_B = csr_max_row_wght(p->spaH);
  ctx.cc_branch = (p->spaH->rows > 0) ? fmax(1.0, (double)p->spaH->p[p->spaH->rows] / p->spaH->rows - 1.0) : 1.0;
  ctx.cc_R[0] = 1;
  for (int l = 1; l <= CC_SPLIT_MAX_LEVELS; l++) {
    const int prev = ctx.cc_R[l - 1];
    ctx.cc_R[l] = (prev > 0 && ctx.cc_B > 0 && prev <= INT_MAX / ctx.cc_B) ? prev * ctx.cc_B : 0;
  }
  ctx.round_R = 1;
 
  if ((p->method & 1) && p->ksub > 0 && ctx.rw_threads > 0) {
    double t_ker0 = get_time_sec();
    ctx.N_global = mzd_nullspace(p->spaH);
    ctx.nu = ctx.N_global ? ctx.N_global->nrows : 0;
    if (p->spaH->rows < ctx.nu || ctx.nu <= 0) {
      if (p->debug & DBG_SUMMARY) {
        fprintf(stderr,
                "# Warning: ksub=%d requested, but m=%d < nu=%d; "
                "falling back to full-matrix RW (ksub=0)\n",
                p->ksub, p->spaH->rows, ctx.nu);
      }
      safe_mzd_free(ctx.N_global);
      ctx.N_global = NULL;
      ctx.nu = 0;
      p->ksub = 0;
      if (!p->nothrottle) { /* full-matrix RW workers need more memory: re-apply the RW thread limit */
        const int rw_cap = rw_thread_cap(p);
        if (p->method == 1 && num_threads > rw_cap) {
          num_threads = rw_cap;
          ctx.num_threads = rw_cap;
        }
        if (ctx.rw_threads > rw_cap) ctx.rw_threads = rw_cap;
        if (p->debug & DBG_TIMING) {
          fprintf(stderr, "# Note: %d threads, %d of them for RW (full-matrix RW)\n", num_threads, ctx.rw_threads);
        }
      }
    } else if (p->debug & DBG_TIMING) {
      fprintf(stderr,
              "# computed ker(H) in %.4fs: nu=%d, n=%d, ksub=%d (eff=%d)\n",
              get_time_sec() - t_ker0, ctx.nu, p->spaH->cols, p->ksub,
              minint(p->ksub, ctx.nu));
    }
    if (ctx.nu > 0) atomic_store(&ctx.rw_rank, p->spaH->cols - ctx.nu); /* rank of H */
  }
 
  /* method 2 starts a few workers only, and the others once a CC round takes longer (see CC_SPAWN_INITIAL) */
  ctx.spawn_initial = (p->method == 2 && !p->nothrottle && num_threads > CC_SPAWN_INITIAL) ? CC_SPAWN_INITIAL
                                                                                           : num_threads;
  if (ctx.spawn_initial < num_threads && (p->debug & DBG_TIMING)) {
    fprintf(stderr, "# Note: CC starts %d of %d threads, the others when a CC round takes longer\n",
            ctx.spawn_initial, num_threads);
  }
 
  if (p->debug & DBG_PROGRESS) print_run_plan(&ctx, init_dmax);
 
  /* Allocate and launch worker threads */
  ctx.threads = malloc(num_threads * sizeof(pthread_t));
  worker_arg_t *args = malloc(num_threads * sizeof(worker_arg_t));
  if (!ctx.threads || !args) ERROR("memory allocation");
  ctx.wargs = args;
 
  for (int i = 0; i < num_threads; i++) {
    args[i].ctx = &ctx;
    args[i].tid = i;
    args[i].irr_ws = NULL; /* (allocated by the worker, see worker_thread_func()) */
    for (int k = 0; k < MAX_W; k++) {
      args[i].min_swei[k] = p->spaH->rows + 1;
    }
  }
  ctx_spawn(&ctx, ctx.spawn_initial);
 
  if (p->method == 1) {
    run_method1_coordinator(&ctx);
  } else if (p->method == 2) {
    run_method2_coordinator(&ctx);
  } else if (p->method == 3) {
    run_method3_coordinator(&ctx);
  } else {
    ERROR("invalid method %d\n", p->method);
  }
 
  /* Signal stop and wait for all workers */
  atomic_store(&ctx.stop_flag, true);
  for (int i = 0; i < ctx.spawned; i++) {
    pthread_join(ctx.threads[i], NULL);
  }
 
  int final_dmin = atomic_load(&ctx.dmin);
  int final_dmax = atomic_load(&ctx.dmax);
  int cc_found = atomic_load(&ctx.cc_found_weight);
  /* A CC codeword gives the exact distance only if found in a round where all lower weights had
   * been analyzed (not with noscan=1 unless dmin=wmax is supplied); see `cc_exact`. */
  const bool cc_exact = (cc_found > 0) && atomic_load(&ctx.cc_exact) &&
                        (final_dmax == 0 || cc_found <= final_dmax);
 
  if (cc_exact) {
    final_dmin = cc_found;
    final_dmax = cc_found;
  } else if (final_dmax > 0 && final_dmin >= final_dmax) {
    final_dmin = final_dmax;
  }
 
  /* why the run ended (e.g., RW ends the run on wmin, unless codewords are collected with outC or maxC) */
  if (p->debug & DBG_SUMMARY) print_stop_reason(&ctx, final_dmax, cc_exact);
 
  /* Confinement profile output (if smax > 0 and CC was run) */
  if (p->smax && p->method >= 2) {
    int max_w_analyzed = (final_dmin > 1) ? (final_dmin - 1) : ((p->wmax > 0) ? p->wmax : 0);
    if (cc_exact) max_w_analyzed = cc_found;
    if (max_w_analyzed > 0) {
      int global_swei[MAX_W];
      for (int i = 0; i < MAX_W; i++) global_swei[i] = p->spaH->rows + 1;
      for (int t = 0; t < ctx.spawned; t++) {
        for (int i = 1; i <= max_w_analyzed; i++) {
          if (args[t].min_swei[i] < global_swei[i]) {
            global_swei[i] = args[t].min_swei[i];
          }
        }
      }
      int skipped = 0;
      if (p->debug & DBG_SUMMARY) {
        for (int i = 1; i <= max_w_analyzed; i++) {
          if (global_swei[i] <= p->spaH->rows) {
            fprintf(stderr, "# w=%d min non-zero syndrome weight %d\n", i, global_swei[i]);
          } else {
            skipped = 1;
          }
        }
      } else {
        fprintf(stderr, "# confinement: ");
        for (int i = 1; i <= max_w_analyzed; i++) {
          if (global_swei[i] <= p->spaH->rows) {
            fprintf(stderr, "%d%s", global_swei[i], i < max_w_analyzed ? "," : "");
          } else {
            skipped = 1;
            fprintf(stderr, "?%s", i < max_w_analyzed ? "," : "");
          }
        }
        fprintf(stderr, "\n");
      }
      if (skipped) {
        fprintf(stderr,
                "# Note: Some weights were skipped in confinement profile. "
                "Try increasing smax (current: %d)\n", p->smax);
      }
    }
  }
 
  long reported_rw_steps = 0;
  if (p->method != 2 && !cc_exact) {
    reported_rw_steps = atomic_load(&ctx.rw_steps_completed);
  }
 
  if (p->debug & DBG_SUMMARY) {
    print_codeword_stats(stderr, p);
    const long long red_cc = atomic_load(&ctx.cw_red_cc), red_rw = atomic_load(&ctx.cw_red_rw);
    if (red_cc > 0 || red_rw > 0) /* (see cw_irreducible(); a codeword found several times is counted each time) */
      fprintf(stderr, "# reducible codewords skipped (irred=1, a codeword plus a disjoint nonzero vector with zero "
              "syndrome): %lld times in CC, %lld times in RW\n", red_cc, red_rw);
    /* RW data for the information-set estimate (see print_rw_infoset_estimate(); also used by the Python wrapper) */
    const long rw_done = atomic_load(&ctx.rw_steps_completed);
    const long rw_unif = atomic_load(&ctx.rw_unif_steps);
    const int rw_rank = atomic_load(&ctx.rw_rank);
    const long rw_timed = atomic_load(&ctx.rw_timed_steps);
    const double t_step = (rw_timed > 0) ? 1e-9 * (double)atomic_load(&ctx.rw_busy_ns) / (double)rw_timed : 0.0;
    const int rw_pairs_l = atomic_load(&ctx.rw_pairs_l);
    if (rw_done > 0 && rw_rank > 0) {
      fprintf(stderr, "# RW information sets: n=%d, rank(H)=%d, steps=%ld (uniform permutations: %ld), ksub=%d, "
              "pairs=%d, %.3g s/step per thread, %d RW threads\n", p->nvar, rw_rank, rw_done, rw_unif, p->ksub,
              rw_pairs_l, t_step, ctx.rw_threads);
      if (rw_pairs_l > 0)
        fprintf(stderr, "# RW pairs (p=2, buckets of l=%d random pivot rows): %.3g pairs checked per step; codewords "
                "recorded with weight up to the upper bound: %lld from pairs, %lld from single columns\n", rw_pairs_l,
                (double)atomic_load(&ctx.rw_pairs_checked) / (double)rw_done, (long long)atomic_load(&ctx.rw_hits2),
                (long long)atomic_load(&ctx.rw_hits1));
      if (rw_unif > 0 && rw_unif < rw_done) /* the codewords found by the uniform and by the window steps */
        print_rw_step_stats(stderr, p);
    }
    /* the RW upper bound is not certified: warn if the hit counts suggest that some codewords are much harder to
     * find than others of the same weight (one pass over the codewords in hash), and then also print the
     * information-set estimate of the probability that RW missed a lighter codeword (not after the stop target
     * dstop was reached: codewords of weight >= dstop are not needed) */
    if (reported_rw_steps > 0 && final_dmin < final_dmax && atomic_load(&ctx.stop_reason) != STOP_DSTOP &&
        codeword_hit_check(stderr, p) > 0 && rw_unif > 0 && rw_rank > 0) {
      print_rw_infoset_estimate(stderr, p->nvar, rw_rank, (final_dmin > 1) ? final_dmin : 1, final_dmax - 1,
                                rw_unif, rw_done, (ctx.rw_threads > 0) ? t_step / ctx.rw_threads : 0.0, rw_pairs_l);
    }
  }
 
  /* Output to stdout: dmin dmax rw_steps */
  printf("%d %d %ld\n", final_dmin, final_dmax, reported_rw_steps);
  fflush(stdout);
 
  /* Codeword export */
  if (p->outC) {
    char comment[256];
    sprintf(comment, "generated by dist_m4ri");
    const long long num_written = nzlist_write(p->outC, comment, p);
    if (p->debug & DBG_SUMMARY) {
      fprintf(stderr, "# exported %lld codewords to %s\n", num_written, p->outC);
    }
  }
 
  if (p->debug & DBG_CWDUMP) {
    cw_vec_t *cw;
    fprintf(stderr, "# all %lld codewords in hash (1-based columns, cnt: number of hits, of which u in RW steps with "
            "uniform permutations and w in RW steps with windows):\n", p->num_cws);
    for (cw = p->codewords; cw != NULL; cw = (cw_vec_t *)(cw->hh.next)) {
      fprintf(stderr, "# cw: [ ");
      for (int i = 0; i < cw->weight; i++) fprintf(stderr, "%d ", 1 + col_orig_index(p, cw->arr[i]));
      fprintf(stderr, "] cnt=%d (u=%d, w=%d)\n", cw->cnt, cw->cnt_u, cw->cnt_w);
    }
  }
 
  /* Cleanup */
  free(rw_lhash);
  safe_mzd_free(ctx.N_global);
  cc_hash_free(cc_hash);
  csr_free(ctx.mHT_cc);
  free(ctx.threads);
  free(args);
  pthread_rwlock_destroy(&ctx.basis_rwlock);
  pthread_mutex_destroy(&ctx.cw_mutex);
 
  var_kill(p);
 
  return 0;
}