Merge branch 'fix/numerics-critical'

Fix the ten defects found by the 2026-09-09 floating-point audit: four
critical, three high, three medium.

Closes #55, #56, #57, #58, #59, #60, #61, #62, #63, #64

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ

src/acquisition.rs

@@ -47,7 +47,38 @@ fn kl_divergence(q: Gaussian, p: Gaussian) -> f64 {
     }
 
     let mean_gap = q.mu() - p.mu();
-    0.5 * (libm::log(var_p / var_q) + (var_q + mean_gap * mean_gap) / var_p - 1.0)
+
+    // Algebraically `0.5 * (ln(var_p/var_q) + (var_q + gap^2)/var_p - 1)`, but
+    // written so that neither term can go negative.
+    //
+    // The direct form cancels against its `- 1.0` for two near-identical
+    // distributions and returns a *negative* divergence — measured, 762 082 of
+    // 3 000 000 near-identical pairs, worst `-5.55e-17`, which is exactly one
+    // ULP of the 1.0. It also loses the answer entirely where it is small:
+    // at `var_q/var_p - 1 = 1e-9` the direct form gives `0.0` where the true
+    // value is `2.5e-19`.
+    //
+    // With `u = var_q/var_p - 1` the variance part is `0.5 * (u - ln(1+u))`,
+    // which is non-negative for every `u > -1`, and the mean part is a square
+    // over a positive variance. Non-negativity is then structural rather than
+    // incidental.
+    let u = var_q / var_p - 1.0;
+    0.5 * u_minus_ln1p(u) + mean_gap * mean_gap / (2.0 * var_p)
+}
+
+/// `u - ln(1 + u)`, without the cancellation that spelling invites.
+///
+/// Both terms are approximately `u` for small `u`, so the subtraction loses
+/// everything just where the result matters. The Taylor series
+/// `u^2/2 - u^3/3 + u^4/4 - ...` is exact in that regime and manifestly
+/// non-negative, since `u^2/2` dominates.
+fn u_minus_ln1p(u: f64) -> f64 {
+    if u.abs() < 1e-4 {
+        let u2 = u * u;
+        u2 * (0.5 - u / 3.0 + u2 / 4.0)
+    } else {
+        u - libm::log1p(u)
+    }
 }
 
 /// Expected information gain of a hypothetical matchup, in nats.
@@ -146,7 +177,7 @@ pub fn expected_information_gain<T: Time, D: Drift<T>>(
 
     let mut gain = 0.0;
 
-    for (ranks, probability) in predict::outcome_distribution(&performances, &margins) {
+    for (ranks, probability) in predict::outcome_distribution(&performances, &margins)? {
         if probability <= NEGLIGIBLE {
             continue;
         }

src/drift.rs

@@ -21,6 +21,23 @@ pub trait Drift<T: Time>: Copy + Debug + Send + Sync {
 ///
 /// For `Time = i64`: variance added is `(to - from) * gamma^2`.
 /// For `Time = Untimed`: elapsed is always 0, so drift is always 0.
+///
+/// # The sign of `gamma` is not meaningful
+///
+/// `gamma` enters only as `gamma * gamma`, so `ConstantDrift(-0.05)` produces
+/// results **bit identical** to `ConstantDrift(0.05)`. That is the same
+/// sign-absorption `HistoryBuilder::sigma`, `HistoryBuilder::beta`,
+/// `Gaussian::from_ms` and `Rating::new` all reject outright.
+///
+/// It is not rejected here because the field is public and positional, so
+/// there is no constructor to intercept — sealing it would break every
+/// `ConstantDrift(x)` in existence for a case whose *resulting model* is
+/// perfectly valid, just not the one a caller writing a minus sign expected.
+///
+/// A non-finite `gamma` is a different matter and **is** rejected:
+/// `History::converge` validates the drift variance each competitor actually
+/// accumulates, which also covers a custom [`Drift`] implementation, and
+/// reports `InferenceError::InvalidParameter`.
 #[derive(Clone, Copy, Debug)]
 pub struct ConstantDrift(pub f64);
 

src/error.rs

@@ -131,6 +131,28 @@ pub enum InferenceError {
     AlreadyRegistered { key: String },
     /// A prediction was given a team with no members.
     EmptyTeam { team: usize },
+    /// The prediction grid cannot resolve the narrowest feature in the matchup.
+    ///
+    /// `predict_outcome` and `predict_ranking` integrate every team's density
+    /// on one shared grid, whose resolution is set by the narrowest sigma (or a
+    /// narrower draw margin). When the widest and narrowest are far enough
+    /// apart, resolving the narrow one across the wide one's support needs more
+    /// nodes than the grid is allowed to hold.
+    ///
+    /// Reported rather than clamped. Clamping is what this replaced, and it
+    /// returned probabilities greater than one — measured, a `P` of 2.79 and a
+    /// `Prediction::total()` of 5.41 — because the trapezoid rule stops
+    /// resolving a density once the step exceeds roughly 1.7 of its sigma.
+    ///
+    /// `predict_win_probabilities` answers the same matchup through adaptive
+    /// quadrature and is accurate here; use it when only the per-team win
+    /// probabilities are needed.
+    GridTooCoarse {
+        /// Nodes required to resolve the narrowest feature.
+        needed: usize,
+        /// Nodes the grid may hold.
+        max: usize,
+    },
     /// A joint posterior was requested where one cannot be formed exactly.
     JointUnavailable { reason: &'static str },
     /// Fewer than two teams were supplied to a prediction.
@@ -219,6 +241,15 @@ impl fmt::Display for InferenceError {
             Self::EmptyTeam { team } => {
                 write!(f, "team {team} has no members")
             }
+            Self::GridTooCoarse { needed, max } => {
+                write!(
+                    f,
+                    "the prediction grid needs {needed} nodes to resolve the narrowest \
+                     team's density across the widest team's support, but may hold only \
+                     {max}; the sigmas in this matchup are too far apart to integrate on \
+                     one grid. Use predict_win_probabilities, which is accurate here"
+                )
+            }
             Self::JointUnavailable { reason } => {
                 write!(f, "no exact joint posterior is available: {reason}")
             }

src/factor/margin.rs

@@ -81,7 +81,7 @@ fn cavity_log_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 {
     // `hypot`, not `sqrt(a^2 + b^2)`: squaring overflows to infinity above a
     // sigma of ~1.3e154 and flushes to zero below ~1.5e-154, and `Gaussian`'s
     // constructors are public so a caller can reach both.
-    let combined_sigma = cavity.sigma().hypot(sigma);
+    let combined_sigma = libm::hypot(cavity.sigma(), sigma);
     let value = ln_pdf(m_obs, cavity.mu(), combined_sigma);
 
     // A degenerate cavity (infinite sigma) is the only way to reach a
@@ -89,7 +89,7 @@ fn cavity_log_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 {
     if value.is_finite() {
         value
     } else {
-        f64::MIN_POSITIVE.ln()
+        libm::log(f64::MIN_POSITIVE)
     }
 }
 

src/factor/trunc.rs

@@ -95,7 +95,7 @@ fn cavity_log_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 {
     if value.is_finite() {
         value
     } else {
-        f64::MIN_POSITIVE.ln()
+        libm::log(f64::MIN_POSITIVE)
     }
 }
 
@@ -203,7 +203,7 @@ mod tests {
             let approx = -0.5 * z * z - z.ln() - (2.0 * std::f64::consts::PI).sqrt().ln();
 
             assert!(
-                got < f64::MIN_POSITIVE.ln(),
+                got < libm::log(f64::MIN_POSITIVE),
                 "mu={mu}: {got} is still stuck on the old clamp floor"
             );
             assert!(

src/gaussian.rs

@@ -18,8 +18,44 @@ pub struct Gaussian {
 
 impl Gaussian {
     /// Construct from mean and standard deviation.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `sigma` is negative. NaN is deliberately allowed through: a
+    /// broken fit produces one, and `converge` reports that as
+    /// `NonFiniteResult` rather than panicking mid-inference.
+    ///
+    /// A negative sigma used to be accepted and returned results **bit
+    /// identical** to its absolute value, because sigma only ever enters as
+    /// `sigma * sigma`. The sign was not rejected and not honoured; it simply
+    /// vanished. That is the same defect `HistoryBuilder::sigma`,
+    /// `HistoryBuilder::beta` and `Member::with_drift_scale` already reject.
+    ///
+    /// # Very small sigma
+    ///
+    /// `pi = 1 / sigma^2` leaves `f64`'s range below about `1.5e-154`, and
+    /// `tau = mu * pi` overflows sooner still — at a threshold that depends on
+    /// `mu`, so there is a band where `pi` is finite and only `tau` is not.
+    /// Both land on the same point-mass representation the `sigma == 0.0`
+    /// branch produces, and a point mass with a non-zero mean has `mu() = NaN`,
+    /// because `tau / pi` is `inf / inf`.
+    ///
+    /// This is not rejected, because `approx` legitimately produces a very
+    /// small truncated sigma and inference must not panic. It is worth knowing
+    /// that such a `Gaussian` is not equal to itself, so two identical
+    /// declarations of one can be reported as conflicting.
     #[must_use]
     pub const fn from_ms(mu: f64, sigma: f64) -> Self {
+        // NaN is admitted on purpose. A broken fit legitimately produces a NaN
+        // sigma — `sqrt` of a negative truncated variance — and the design is
+        // to propagate that to `converge`'s `NonFiniteResult` guard, not to
+        // panic inside inference. Rejecting it here turned that reporting path
+        // into a crash, which two tests caught immediately.
+        assert!(
+            sigma >= 0.0 || sigma.is_nan(),
+            "sigma must not be negative; it is only ever squared, so a negative \
+             value would silently behave as its absolute value"
+        );
         if sigma == f64::INFINITY {
             Self { pi: 0.0, tau: 0.0 }
         } else if sigma == 0.0 {
@@ -120,7 +156,25 @@ impl Gaussian {
         }
     }
 
+    /// How far this Gaussian moved from `other`, as `(|d mu|, |d sigma|)`.
+    ///
+    /// Identical messages have not moved, whatever their parameters, and that
+    /// case is answered in natural space before touching `mu()`/`sigma()`. An
+    /// improper message has `pi == 0`, so `sigma()` is infinite — and
+    /// `inf - inf` is NaN, a NaN *change* for a message that did not change at
+    /// all. (`mu()` is guarded and returns 0.0 here, so the mean component was
+    /// never the problem; the sigma component alone produced `(0.0, NaN)`.)
+    ///
+    /// That is reachable in ordinary inference: once a pairing is more than
+    /// about nine cavity-sigma apart the truncation is a no-op, `trunc / cavity`
+    /// is exactly the identity message, and the chain compares one identity
+    /// against another. Before this guard that produced `(0.0, NaN)`, which
+    /// silently disabled the sigma half of the convergence test.
     pub(crate) fn delta(&self, other: Gaussian) -> (f64, f64) {
+        if self.pi == other.pi && self.tau == other.tau {
+            return (0.0, 0.0);
+        }
+
         (
             (self.mu() - other.mu()).abs(),
             (self.sigma() - other.sigma()).abs(),
@@ -256,6 +310,42 @@ impl ops::Div<Gaussian> for Gaussian {
 
 #[cfg(test)]
 mod tests {
+    /// A message that did not change must report no change, even when it is
+    /// improper. `mu()` of an improper Gaussian is `0/0 = NaN` and `sigma()` is
+    /// infinite, so the mean/sigma form reported `(NaN, NaN)` for two identical
+    /// identity messages — which silently disabled the sigma half of the
+    /// convergence test in `run_chain`.
+    #[test]
+    fn delta_of_two_identical_improper_messages_is_zero() {
+        let improper = crate::N_INF;
+        // `mu()` is guarded and returns 0.0 for an improper Gaussian, so the
+        // mean component was always fine. The NaN came from the sigma
+        // component alone: `inf - inf`. The pre-fix value was `(0.0, NaN)`.
+        assert!(improper.sigma().is_infinite(), "premise: sigma is infinite");
+        assert_eq!(improper.mu(), 0.0, "premise: mu is guarded, not NaN");
+        assert!(
+            (improper.sigma() - improper.sigma()).is_nan(),
+            "premise: the unguarded sigma difference is NaN"
+        );
+        assert_eq!(improper.delta(improper), (0.0, 0.0));
+    }
+
+    #[test]
+    fn delta_of_identical_proper_messages_is_zero() {
+        let g = Gaussian::from_ms(25.0, 8.0);
+        assert_eq!(g.delta(g), (0.0, 0.0));
+    }
+
+    /// The shortcut must not swallow a real difference.
+    #[test]
+    fn delta_still_measures_a_real_move() {
+        let a = Gaussian::from_ms(25.0, 8.0);
+        let b = Gaussian::from_ms(26.0, 9.0);
+        let (dmu, dsigma) = a.delta(b);
+        assert!((dmu - 1.0).abs() < 1e-12, "{dmu}");
+        assert!((dsigma - 1.0).abs() < 1e-12, "{dsigma}");
+    }
+
     use super::*;
 
     #[test]

src/history.rs

@@ -1,4 +1,9 @@
-use std::{borrow::Borrow, collections::HashMap, hash::Hash, marker::PhantomData};
+use std::{
+    borrow::Borrow,
+    collections::{BTreeMap, HashMap},
+    hash::Hash,
+    marker::PhantomData,
+};
 
 use crate::{
     BETA, GAMMA, Index, MU, P_DRAW, SIGMA,
@@ -261,7 +266,15 @@ struct ResolvedTerms {
     contrast: Vec<f64>,
     /// Coefficients of competitors the slice has never seen, keyed by their
     /// rendering. Independent of everything in the slice by construction.
-    unseen: HashMap<String, f64>,
+    ///
+    /// A `BTreeMap` rather than a `HashMap`, and that is load-bearing. These
+    /// coefficients are summed, addition is not associative, and Rust seeds its
+    /// default hasher per process — so iterating a `HashMap` here made
+    /// `posterior_of` return different bits run to run on identical input.
+    /// Measured over 40 processes: two distinct sigma bit patterns, and five
+    /// distinct values from `expected_variance_reduction` spanning ~7 ULP.
+    /// Ordered iteration makes the sum reproducible.
+    unseen: BTreeMap<String, f64>,
     mean: f64,
 }
 
@@ -947,6 +960,40 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     ///
     /// Returns the matrix, each competitor's row at its latest appearance, and
     /// the row at each `(competitor, slice)` for time-addressed queries.
+    /// Drift below which two consecutive appearances are one latent variable.
+    ///
+    /// The rule used to be `drift <= 0.0` exactly, and everything above it got an
+    /// explicit `1.0 / drift` precision. That is a representation the matrix cannot
+    /// hold: at `drift = 1e-16` the entry is `1e16`, and `1e16 + 0.28` rounds back
+    /// to `1e16`, so the prior and the event contrasts are annihilated in the
+    /// stored `f64` before the factorisation ever runs. Measured, `drift_scale =
+    /// 1e-10` returned a posterior variance **12 000x too small** — a 111x
+    /// overconfident interval — as `Ok`, and the band just above it returned a
+    /// misleading `JointUnavailable`.
+    ///
+    /// Solved exactly in high precision the same system is perfectly well
+    /// conditioned: it converges smoothly onto the collapsed value and is flat from
+    /// `1e-16` down to `1e-40`. So this is a representation problem, not a
+    /// conditioning one, and scaling cannot fix it — symmetric (Jacobi)
+    /// equilibration was measured **30x worse**, because the information is already
+    /// gone from the assembled matrix by the time a solver sees it.
+    ///
+    /// The threshold balances the two errors that trade off here. Ignoring a real
+    /// drift costs roughly `drift / V`; representing one costs roughly
+    /// `EPSILON * V / drift`, since `1 / drift` swamps the other precisions in the
+    /// row. They cross at `drift ~ V * sqrt(EPSILON)`, which is what this returns.
+    /// `V` is the competitor's own prior variance, so the threshold follows the
+    /// scale each competitor is actually measured on.
+    ///
+    /// Ordinary drift is far above this and is unaffected: the crate's default
+    /// `gamma = 25/300` accumulates `0.0069` per unit time against a threshold of
+    /// `1.0e-6` at the default prior.
+    fn collapse_threshold(prior_variance: f64) -> f64 {
+        // sqrt(f64::EPSILON), as a const rather than a runtime sqrt.
+        const SQRT_EPSILON: f64 = 1.490_116_119_384_765_6e-8;
+        prior_variance * SQRT_EPSILON
+    }
+
     fn time_expanded_joint(&self) -> TimeExpanded {
         let mut latest: HashMap<Index, (usize, usize)> = HashMap::new();
         let mut at_slice: HashMap<(Index, usize), usize> = HashMap::new();
@@ -969,8 +1016,9 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
                     }
                     Some(&prev) => {
                         let drift = rating.drift_variance_for_elapsed(elapsed);
-                        if drift <= 0.0 {
-                            // No drift: the same latent skill, not two.
+                        if drift <= Self::collapse_threshold(rating.prior.variance()) {
+                            // No drift, or too little to represent: the same
+                            // latent skill, not two. See `collapse_threshold`.
                             prev
                         } else {
                             let row = n;
@@ -1032,7 +1080,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         K: std::fmt::Debug,
     {
         let mut contrast = vec![0.0; width];
-        let mut unseen: HashMap<String, f64> = HashMap::new();
+        let mut unseen: BTreeMap<String, f64> = BTreeMap::new();
         let mut mean = 0.0;
 
         for (member, (key, coefficient)) in terms.iter().enumerate() {
@@ -1278,8 +1326,10 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
 
         let cholesky = crate::joint::Cholesky::factor(lambda, width).ok_or(
             InferenceError::JointUnavailable {
-                reason: "the precision matrix is not positive-definite, which means \
-                         a competitor has neither a proper prior nor any evidence",
+                reason: "the precision matrix is not positive-definite; the usual \
+                         cause is a competitor with neither a proper prior nor any \
+                         evidence, but an extreme prior or drift can also make the \
+                         assembled matrix indefinite in floating point",
             },
         )?;
 
@@ -1502,7 +1552,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         Ok(Prediction::new(crate::predict::outcome_distribution(
             &performances,
             &self.margins(&sizes),
-        )))
+        )?))
     }
 
     /// Probability of one specific finishing order.
@@ -1542,11 +1592,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         }
 
         let (performances, sizes) = self.performances(teams)?;
-        Ok(crate::predict::ranking_probability(
-            &performances,
-            &self.margins(&sizes),
-            ranks,
-        ))
+        crate::predict::ranking_probability(&performances, &self.margins(&sizes), ranks)
     }
 
     /// Run the full forward+backward convergence loop to a fixed point.
@@ -1611,6 +1657,30 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
 
         let opts = self.convergence;
 
+        // Drift is the one model parameter with no boundary check, because
+        // `HistoryBuilder::drift` is generic over `Drift<T>` and cannot inspect
+        // an arbitrary implementation. Validate what it actually produces
+        // instead, which also covers a custom impl.
+        //
+        // `ConstantDrift` returns `elapsed * gamma * gamma`, so a negative
+        // gamma is squared away: measured, `ConstantDrift(-0.0833)` gave
+        // results **bit identical** to `+0.0833`, the same sign-absorption
+        // defect already rejected for `sigma` and `beta`. A non-finite gamma
+        // poisons every posterior derived from it.
+        for slice in &self.time_slices {
+            for (agent, elapsed) in slice.appearances() {
+                let drift = self.agents[agent]
+                    .rating
+                    .drift_variance_for_elapsed(elapsed);
+                if !drift.is_finite() || drift < 0.0 {
+                    return Err(InferenceError::InvalidParameter {
+                        name: "drift variance",
+                        value: drift,
+                    });
+                }
+            }
+        }
+
         if self.time_slices.is_empty() {
             return Ok(ConvergenceReport {
                 iterations: 0,
@@ -1801,7 +1871,16 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         // batched and succeeded, order-dependently, when fed one at a time.
         // Checked before anything mutates, so a rejected batch leaves the
         // history untouched.
-        for (agent, batch) in &priors {
+        // Sorted, not `HashMap` order. This loop returns on the FIRST conflict
+        // it finds, so hash order decided *which* competitor the error blamed:
+        // measured, 15 different competitors named across 40 runs on identical
+        // input. The error fired every time — only its content was a lottery,
+        // which makes it unreproducible and sends a reader after the wrong key.
+        let mut conflict_scan: Vec<Index> = priors.keys().copied().collect();
+        conflict_scan.sort_unstable();
+
+        for agent in &conflict_scan {
+            let batch = priors[agent];
             let held = self.declared.get(agent).copied().unwrap_or_default();
 
             if let (Some(existing), Some(new)) = (held.prior, batch.prior) {

src/lib.rs

@@ -218,6 +218,12 @@ const HALF_LINE_WINDOW: f64 = 10.0;
 /// four-term series is good to ~1e-10 by here, so the two are at their closest
 /// agreement around this point. Below it the subtraction is exact; above it the
 /// series is.
+/// `alpha / width` past which the tie branch's `v^2 - u` has lost too many
+/// digits to trust, and the narrow-window form takes over.
+///
+/// The subtraction retains about `(width / alpha)^2 / EPSILON` of its
+/// precision, so this is the ratio at which that falls below roughly 1e-6.
+const NARROW_WINDOW_RATIO: f64 = 2.0e4;
 const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0;
 
 pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0);
@@ -476,10 +482,72 @@ fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
 fn half_line_truncation(alpha: f64) -> (f64, f64) {
     let inv = alpha.recip();
     let inv_sq = inv * inv;
-    let gap = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
+    let b = 2.0 - inv_sq * (10.0 - 74.0 * inv_sq);
+    let gap = inv * (1.0 - inv_sq * b);
     let v = alpha + gap;
 
-    (v, v * gap)
+    // Returns `1 - w`, not `w`, and that is the whole point of this shape.
+    //
+    // `w` tends to 1 out here, so a caller forming `1 - w` loses about
+    // `log10(alpha^2)` digits: measured against the exact truncated variance,
+    // `1 - w` came back with 8.9e-5 relative error at alpha = 1e6 and **0.0**
+    // from alpha = 1e8 — where the true value is 1e-16 and perfectly
+    // representable. `sigma * (1 - w).sqrt()` was then exactly zero, and
+    // `from_ms(mu, 0.0)` is a point mass whose `mu()` is `inf/inf = NaN`.
+    //
+    // Expanding `1 - v*gap` symbolically removes the subtraction: with
+    // `alpha*gap = 1 - inv^2*b`, the leading ones cancel on paper instead of in
+    // floating point, leaving `inv^2` times a bracket that tends to 1. Measured
+    // exact — 0.0 relative error — from alpha = 1e3 to 1e8.
+    let one_minus_w = inv_sq
+        * ((1.0 - inv_sq * (10.0 - 74.0 * inv_sq)) + 2.0 * inv_sq * b - inv_sq * inv_sq * b * b);
+
+    (v, one_minus_w)
+}
+
+/// Truncation to a *narrow* window `[alpha, alpha + d]`, as `(v, 1 - w)`.
+///
+/// The tie branch forms `w` from `v^2 - u`, and both grow as `alpha^2` while
+/// their difference stays `O(1)`. Far enough into the tail that subtraction has
+/// nothing left: measured at `alpha = 1e6` with a window of `1e-6` it kept four
+/// significant digits and returned `1 - w = -2.4e-4` where the truth is
+/// `+2.8e-13`, so `sqrt` of it was NaN. One step earlier it was quietly wrong
+/// instead — `1 - w = 1.0` exactly, a truncation reported as a no-op, where the
+/// truth was `5e-17`.
+///
+/// The existing half-line escape hatch does not cover it, because that keys on
+/// `alpha * d >= HALF_LINE_WINDOW` — how many window-widths from the mean the
+/// window sits — and a *narrow* window fails that however deep it is.
+///
+/// Over a narrow window the density is `exp(-t*s - s^2 d^2 / 2)` in
+/// `x = alpha + s*d`, with `t = alpha * d`. Dropping the `d^2` term leaves a
+/// truncated exponential on `[0, 1]`, whose mean and variance are closed forms.
+/// So `v = alpha + d*m(t)` and `1 - w = d^2 * V(t)`, with no subtraction of
+/// large quantities anywhere.
+///
+/// Measured against high-precision quadrature over `alpha` in `[1e2, 1e9]`:
+/// `v` exact to 4e-10 or better, `1 - w` to 4e-10 across the region this is
+/// used in.
+fn narrow_window_truncation(alpha: f64, d: f64) -> (f64, f64) {
+    let t = alpha * d;
+
+    // `m` and `V` are the mean and variance of a truncated exponential on
+    // [0, 1] with rate `t`, both of which cancel as `t -> 0`. The series is
+    // their limit (1/2 and 1/12, a uniform window) with the leading correction.
+    let (m, v_s) = if t < 1e-3 {
+        (
+            0.5 - t / 12.0 + t * t * t / 720.0,
+            1.0 / 12.0 - t * t / 240.0,
+        )
+    } else {
+        let em1 = libm::expm1(t);
+        (
+            1.0 / t - 1.0 / em1,
+            1.0 / (t * t) - (em1 + 1.0) / (em1 * em1),
+        )
+    };
+
+    (alpha + d * m, d * d * v_s)
 }
 
 fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
@@ -507,7 +575,7 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
             (v, v - alpha)
         };
 
-        (v, v * gap)
+        (v, 1.0 - v * gap)
     } else {
         // v is odd in mu and w is even, so fold to mu <= 0. Both truncation
         // points then sit in the upper tail, where the scaled form applies.
@@ -523,9 +591,22 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
         // Once the window sits many of its own widths into the tail it is
         // indistinguishable from a half-line, so the asymptotic covers it with
         // no subtraction at all.
-        if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * (beta - alpha) >= HALF_LINE_WINDOW {
-            let (v, w) = half_line_truncation(alpha);
-            return (if flipped { -v } else { v }, w);
+        let width = beta - alpha;
+
+        if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * width >= HALF_LINE_WINDOW {
+            let (v, one_minus_w) = half_line_truncation(alpha);
+            return (if flipped { -v } else { v }, one_minus_w);
+        }
+
+        // A narrow window deep in the tail: too narrow for the half-line above,
+        // too deep for the subtraction below. The direct form keeps roughly
+        // `1 / (alpha/width)^2` of its digits, so the crossover is on that
+        // ratio rather than on either quantity alone — and the approximation is
+        // most accurate exactly where the subtraction is worst, since both
+        // improve as the window narrows.
+        if alpha > 0.0 && alpha > NARROW_WINDOW_RATIO * width {
+            let (v, one_minus_w) = narrow_window_truncation(alpha, width);
+            return (if flipped { -v } else { v }, one_minus_w);
         }
 
         let (v, u) = if alpha > 0.0 {
@@ -548,17 +629,23 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
             )
         };
 
-        let w = -(u - v.powi(2));
+        // `1 - w` where `w = v^2 - u`. Both `v^2` and `u` grow as alpha^2 while
+        // their difference stays O(1), so this subtraction is the one place the
+        // tie branch can still lose everything — see the escape hatch above,
+        // which is what keeps the far tail away from it.
+        let one_minus_w = 1.0 + u - v.powi(2);
 
-        (if flipped { -v } else { v }, w)
+        (if flipped { -v } else { v }, one_minus_w)
     }
 }
 
 fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
-    let (v, w) = v_w(mu, sigma, margin, tie);
+    // `v_w` returns `1 - w` rather than `w`: forming the difference here is
+    // what destroyed the truncated variance in the far tail.
+    let (v, one_minus_w) = v_w(mu, sigma, margin, tie);
 
     let mu_trunc = mu + sigma * v;
-    let sigma_trunc = sigma * (1.0 - w).sqrt();
+    let sigma_trunc = sigma * one_minus_w.sqrt();
 
     (mu_trunc, sigma_trunc)
 }
@@ -569,13 +656,34 @@ pub(crate) fn approx(n: Gaussian, margin: f64, tie: bool) -> Gaussian {
     Gaussian::from_ms(mu, sigma)
 }
 
+/// Componentwise maximum that **propagates** NaN rather than dropping it.
+///
+/// Every caller folds this as `tuple_max(accumulator, new)`. A plain `>`
+/// comparison is false against NaN, so a NaN accumulator would be replaced by
+/// the next finite delta and the breakdown would vanish — leaving `step_is_finite`
+/// to pass on a fit that is already NaN. Because the fold runs over a `HashMap`,
+/// whether that happened depended on per-process hash order: measured, a NaN fit
+/// was reported as `converged: true` in 16 of 30 runs on identical input.
+///
+/// `f64::max` is not a substitute: it also ignores NaN by design, which is the
+/// same defect wearing a standard-library name.
 pub(crate) fn tuple_max(v1: (f64, f64), v2: (f64, f64)) -> (f64, f64) {
     (
-        if v1.0 > v2.0 { v1.0 } else { v2.0 },
-        if v1.1 > v2.1 { v1.1 } else { v2.1 },
+        max_propagating_nan(v1.0, v2.0),
+        max_propagating_nan(v1.1, v2.1),
     )
 }
 
+fn max_propagating_nan(a: f64, b: f64) -> f64 {
+    if a.is_nan() || b.is_nan() {
+        f64::NAN
+    } else if a > b {
+        a
+    } else {
+        b
+    }
+}
+
 pub(crate) fn tuple_gt(t: (f64, f64), e: f64) -> bool {
     t.0 > e || t.1 > e
 }
@@ -652,6 +760,13 @@ pub(crate) fn sort_time<T: Copy + Ord>(xs: &[T], reverse: bool) -> Vec<usize> {
 /// Panics if fewer than two rating groups are supplied, or if any group is
 /// empty — match quality is a property of a contest between at least two
 /// non-empty sides.
+///
+/// Also panics with "cannot invert a singular matrix" when every rating has
+/// zero sigma *and* `beta` is zero. Nothing is then uncertain, so there is no
+/// distribution to take the quality of; `Gaussian::from_ms(mu, 0.0)` is a point
+/// mass and its `mu()` is not even well defined. Documented rather than
+/// converted, because the input has no meaningful answer rather than an
+/// awkward one.
 #[must_use]
 pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
     assert!(
@@ -717,13 +832,141 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
     let end = &rotated_a_matrix * &mean_matrix;
 
     let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant();
-    let s_arg = ata.determinant() / middle.determinant();
 
-    libm::exp(e_arg) * s_arg.sqrt()
+    // `sqrt(det(ata) / det(middle))`, taken in log space. Both determinants are
+    // products of `k - 1` diagonal entries, so they leave `f64`'s range long
+    // before their ratio does: measured at the crate defaults, 150 groups was
+    // correct at `8.45e-53`, 200 returned `0`, and 250 returned `NaN` where the
+    // true value is `9.51e-88`. With a small beta it is sharper still — at
+    // `sigma = beta = 1e-3`, 60 groups returned `NaN` against a true `1.32e-9`.
+    //
+    // The ratio is what the answer needs and it is representable throughout, so
+    // the intermediates are the only thing that ever overflowed.
+    let ln_s_arg = ata.ln_abs_determinant() - middle.ln_abs_determinant();
+
+    libm::exp(e_arg + 0.5 * ln_s_arg)
 }
 
 #[cfg(test)]
 mod tests {
+    /// The truncated variance must stay a variance across every branch, and
+    /// the branches must agree where they meet.
+    ///
+    /// `v_w` now has three regimes for a tie — half-line, narrow-window, and
+    /// the direct subtraction — and a misplaced crossover between them is the
+    /// failure mode this guards. A jump at a boundary is visible here even
+    /// though the absolute values are not pinned.
+    #[test]
+    fn truncated_variance_is_continuous_across_the_tie_branches() {
+        for &alpha in &[50.0, 99.0, 100.0, 101.0, 1e3, 1e5, 1e6] {
+            // Sweep the window width across NARROW_WINDOW_RATIO and the
+            // half-line threshold, which sit at different widths per alpha.
+            let mut previous: Option<(f64, f64)> = None;
+            let mut width = alpha / (NARROW_WINDOW_RATIO * 100.0);
+            while width < 40.0 / alpha {
+                // mu = 0 puts the window at [-margin, margin]; shift it out to
+                // `alpha` by moving the mean instead.
+                let margin = width * 0.5;
+                let mu = -(alpha + width * 0.5);
+                let (v, one_minus_w) = v_w(mu, 1.0, margin, true);
+
+                assert!(v.is_finite(), "alpha {alpha}, width {width:e}: v = {v}");
+                assert!(
+                    one_minus_w.is_finite() && one_minus_w > 0.0 && one_minus_w <= 1.0,
+                    "alpha {alpha}, width {width:e}: 1 - w = {one_minus_w:e} is not a variance"
+                );
+
+                if let Some((pv, pw)) = previous {
+                    // Consecutive widths differ by 2x, so the moments may not
+                    // differ by more than a small multiple of that.
+                    assert!(
+                        one_minus_w / pw < 32.0 && pw / one_minus_w < 32.0,
+                        "alpha {alpha}: 1 - w jumped from {pw:e} to {one_minus_w:e} \
+                         at width {width:e} — a branch boundary is misplaced"
+                    );
+                    assert!(
+                        (v - pv).abs() <= 8.0 * width.max(1e-12) + 1e-9 * v.abs(),
+                        "alpha {alpha}: v jumped from {pv} to {v} at width {width:e}"
+                    );
+                }
+                previous = Some((v, one_minus_w));
+                width *= 2.0;
+            }
+        }
+    }
+
+    /// The narrow-window form against high-precision quadrature.
+    ///
+    /// These are the inputs where the direct `v^2 - u` subtraction had four
+    /// significant digits left and returned a negative variance.
+    #[test]
+    fn narrow_window_truncation_matches_quadrature() {
+        for &(alpha, d, expect_v, expect_w) in &[
+            (1e6, 2e-6, 1_000_000.000_000_687, 2.759_383_390_335_666e-13),
+            (1e4, 1e-6, 10_000.000_000_499_167, 8.333_291_666_831_727e-14),
+            (
+                1e3,
+                1e-5,
+                1_000.000_004_991_666_6,
+                8.333_291_666_803_818e-12,
+            ),
+        ] {
+            let (v, one_minus_w) = narrow_window_truncation(alpha, d);
+            assert!(
+                ((v - expect_v) / expect_v).abs() < 1e-12,
+                "alpha {alpha:e}: v = {v}, want {expect_v}"
+            );
+            assert!(
+                ((one_minus_w - expect_w) / expect_w).abs() < 1e-8,
+                "alpha {alpha:e}: 1 - w = {one_minus_w:e}, want {expect_w:e}"
+            );
+        }
+    }
+
+    /// A NaN must survive the fold from ANY position, not only the last.
+    ///
+    /// The fold runs over a `HashMap`, so "last" is per-process hash order. The
+    /// end-to-end symptom was a NaN fit reported as `converged: true` in 16 of
+    /// 30 runs on identical input; these three cases are the deterministic form
+    /// of that, so a regression cannot hide behind a lucky seed.
+    #[test]
+    fn tuple_max_propagates_a_nan_from_any_position() {
+        let nan = (f64::NAN, f64::NAN);
+        let small = (1e-9, 1e-9);
+        let big = (1e-3, 1e-3);
+
+        // NaN last.
+        let step = tuple_max(tuple_max(big, small), nan);
+        assert!(!step_is_finite(step), "NaN last: {step:?}");
+
+        // NaN middle.
+        let step = tuple_max(tuple_max(big, nan), small);
+        assert!(!step_is_finite(step), "NaN middle: {step:?}");
+
+        // NaN first — the case a plain `>` comparison drops.
+        let step = tuple_max(tuple_max(nan, big), small);
+        assert!(!step_is_finite(step), "NaN first: {step:?}");
+    }
+
+    /// `f64::max` would pass the test above's first two cases and fail the
+    /// third, so pin that it is not what we use.
+    #[test]
+    fn tuple_max_is_not_f64_max() {
+        assert!(
+            f64::max(f64::NAN, 1.0) == 1.0,
+            "premise: f64::max drops NaN"
+        );
+        let (a, _) = tuple_max((f64::NAN, 0.0), (1.0, 0.0));
+        assert!(a.is_nan(), "tuple_max must not drop what f64::max drops");
+    }
+
+    /// Ordinary values are unaffected.
+    #[test]
+    fn tuple_max_still_takes_the_larger_component() {
+        assert_eq!(tuple_max((1.0, 5.0), (3.0, 2.0)), (3.0, 5.0));
+        assert_eq!(tuple_max((3.0, 2.0), (1.0, 5.0)), (3.0, 5.0));
+    }
+
     use ::approx::assert_ulps_eq;
 
     use super::*;

src/matrix.rs

@@ -91,6 +91,29 @@ impl Lu {
         det
     }
 
+    /// `ln |det|`, accumulated term by term rather than multiplied out.
+    ///
+    /// The determinant of an `n x n` Gram matrix is a product of `n` diagonal
+    /// entries, so it leaves `f64`'s range long before the quantities built
+    /// from it do. `quality()` only ever wants a *ratio* of two determinants,
+    /// and that ratio is perfectly representable while the determinants
+    /// themselves are not — measured, at 250 rating groups both overflow and
+    /// the ratio came back `NaN` where the true answer is `9.51e-88`.
+    ///
+    /// Returns `-inf` for a singular matrix, so `exp` of it is zero.
+    fn ln_abs_determinant(&self) -> f64 {
+        if self.sign == 0.0 {
+            return f64::NEG_INFINITY;
+        }
+
+        let mut acc = 0.0;
+        for i in 0..self.n {
+            acc += libm::log(self.lu[i * self.n + i].abs());
+        }
+
+        acc
+    }
+
     /// Solve `Ax = b` for a single column of the identity, giving one column
     /// of the inverse.
     fn solve_column(&self, col: usize, out: &mut [f64]) {
@@ -157,6 +180,24 @@ impl Matrix {
         Lu::decompose(self).determinant()
     }
 
+    /// `ln |det|` of a square matrix; `-inf` when singular.
+    ///
+    /// See [`Lu::ln_abs_determinant`] for why a ratio of determinants must be
+    /// taken this way.
+    pub fn ln_abs_determinant(&self) -> f64 {
+        assert_eq!(
+            self.width, self.height,
+            "determinant requires a square matrix, got {}x{}",
+            self.height, self.width
+        );
+
+        if self.width == 0 {
+            return 0.0;
+        }
+
+        Lu::decompose(self).ln_abs_determinant()
+    }
+
     /// Matrix inverse via LU decomposition.
     ///
     /// # Panics

src/predict.rs

@@ -21,7 +21,7 @@
 //! would have made every `predict_*` call return a slightly different number,
 //! which is not a property a rating library should have.
 
-use crate::{Gaussian, quadrature};
+use crate::{Gaussian, InferenceError, quadrature};
 
 /// Teams beyond this count make the outcome enumeration impractical.
 ///
@@ -52,6 +52,10 @@ const WIN_TOLERANCE: f64 = 1e-8;
 /// point where refining stops helping.
 const MIN_GRID_POINTS: usize = 8_192;
 const MAX_GRID_POINTS: usize = 262_144;
+/// Nodes requested across the narrowest feature the recursion must resolve.
+const NODES_PER_FEATURE: f64 = 12.0;
+/// Nodes below which the trapezoid rule stops resolving that feature at all.
+const MIN_NODES_PER_FEATURE: f64 = 4.0;
 
 /// How many standard deviations of support the grid and integrals cover.
 ///
@@ -152,7 +156,7 @@ pub(crate) fn win_probabilities(perf: &[Gaussian], margins: &Margins) -> Vec<f64
 /// Resolution is set by the *smallest* feature in play — the narrowest sigma,
 /// or a draw margin narrower still — because that is what the recursion has to
 /// resolve. A grid sized off the widest team would step over the narrow one.
-fn grid_shape(perf: &[Gaussian], margins: &Margins) -> (f64, f64, usize) {
+fn grid_shape(perf: &[Gaussian], margins: &Margins) -> Result<(f64, f64, usize), InferenceError> {
     let lo = perf
         .iter()
         .map(|g| g.mu() - SUPPORT_SIGMAS * g.sigma())
@@ -175,18 +179,36 @@ fn grid_shape(perf: &[Gaussian], margins: &Margins) -> (f64, f64, usize) {
 
     let feature = narrowest.min(smallest_margin);
     let wanted = if feature.is_finite() && feature > 0.0 {
-        ((hi - lo) / (feature / 12.0)).ceil()
+        ((hi - lo) / (feature / NODES_PER_FEATURE)).ceil()
     } else {
         MIN_GRID_POINTS as f64
     };
 
-    let points = if wanted.is_finite() {
-        (wanted as usize).clamp(MIN_GRID_POINTS, MAX_GRID_POINTS)
-    } else {
-        MIN_GRID_POINTS
-    };
+    if !wanted.is_finite() {
+        return Ok((lo, hi, MIN_GRID_POINTS));
+    }
 
-    (lo, hi, points)
+    // Report rather than clamp. Clamping is what this replaced: it silently
+    // handed the recursion a grid too coarse for the narrowest density, and the
+    // trapezoid rule then returned probabilities greater than one — measured, a
+    // `P` of 2.79 and a total of 5.41. Trapezoid error on a Gaussian is
+    // `~exp(-2 pi^2 (sigma/h)^2)`, which is 1e-12 at `h/sigma = 0.86` and O(1)
+    // by `h/sigma = 17`, so the cliff is sharp and there is no useful answer on
+    // the far side of it.
+    //
+    // The floor is `MIN_NODES_PER_FEATURE` rather than the `NODES_PER_FEATURE`
+    // asked for, because the request carries a large margin: measured accurate
+    // to 2.2e-12 at 1.4 nodes per sigma, and wrong by 1.2e-3 at 0.7.
+    let needed = wanted as usize;
+    let floor = ((hi - lo) / (feature / MIN_NODES_PER_FEATURE)).ceil();
+    if floor.is_finite() && floor as usize > MAX_GRID_POINTS {
+        return Err(InferenceError::GridTooCoarse {
+            needed,
+            max: MAX_GRID_POINTS,
+        });
+    }
+
+    Ok((lo, hi, needed.clamp(MIN_GRID_POINTS, MAX_GRID_POINTS)))
 }
 
 /// Densities of each team sampled on the shared grid.
@@ -198,8 +220,8 @@ struct Sampled {
 }
 
 impl Sampled {
-    fn new(perf: &[Gaussian], margins: &Margins) -> Self {
-        let (lo, hi, points) = grid_shape(perf, margins);
+    fn new(perf: &[Gaussian], margins: &Margins) -> Result<Self, InferenceError> {
+        let (lo, hi, points) = grid_shape(perf, margins)?;
         let step = (hi - lo) / (points - 1) as f64;
         let density = perf
             .iter()
@@ -209,12 +231,12 @@ impl Sampled {
                     .collect()
             })
             .collect();
-        Self {
+        Ok(Self {
             lo,
             step,
             points,
             density,
-        }
+        })
     }
 
     fn node(&self, i: usize) -> f64 {
@@ -317,9 +339,12 @@ fn events(n: usize, strict_only: bool) -> Vec<(Vec<usize>, Vec<bool>)> {
 ///
 /// Orders that differ only *within* a tied group describe the same finishing
 /// order, so their probabilities are summed into one entry.
-pub(crate) fn outcome_distribution(perf: &[Gaussian], margins: &Margins) -> Vec<(Vec<u32>, f64)> {
+pub(crate) fn outcome_distribution(
+    perf: &[Gaussian],
+    margins: &Margins,
+) -> Result<Vec<(Vec<u32>, f64)>, InferenceError> {
     let n = perf.len();
-    let sampled = Sampled::new(perf, margins);
+    let sampled = Sampled::new(perf, margins)?;
 
     let mut aggregated: Vec<(Vec<u32>, f64)> = Vec::new();
     for (order, tied) in events(n, margins.all_zero()) {
@@ -332,7 +357,7 @@ pub(crate) fn outcome_distribution(perf: &[Gaussian], margins: &Margins) -> Vec<
     }
 
     aggregated.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
-    aggregated
+    Ok(aggregated)
 }
 
 /// All permutations of `items`.
@@ -397,9 +422,13 @@ fn orders_for_groups(groups: &[Vec<usize>]) -> Vec<(Vec<usize>, Vec<bool>)> {
 /// Ties in `ranks` mean the tied teams may finish in any internal order, so
 /// this sums the orders consistent with the requested ranking rather than
 /// picking one.
-pub(crate) fn ranking_probability(perf: &[Gaussian], margins: &Margins, ranks: &[u32]) -> f64 {
+pub(crate) fn ranking_probability(
+    perf: &[Gaussian],
+    margins: &Margins,
+    ranks: &[u32],
+) -> Result<f64, InferenceError> {
     let n = perf.len();
-    let sampled = Sampled::new(perf, margins);
+    let sampled = Sampled::new(perf, margins)?;
 
     let mut distinct: Vec<u32> = ranks.to_vec();
     distinct.sort_unstable();
@@ -410,10 +439,10 @@ pub(crate) fn ranking_probability(perf: &[Gaussian], margins: &Margins, ranks: &
         .map(|&r| (0..n).filter(|&i| ranks[i] == r).collect())
         .collect();
 
-    orders_for_groups(&groups)
+    Ok(orders_for_groups(&groups)
         .iter()
         .map(|(order, tied)| order_probability(margins, &sampled, order, tied))
-        .sum()
+        .sum())
 }
 
 /// A distribution over the ways a contest could finish.
@@ -604,7 +633,7 @@ mod tests {
             ),
         ] {
             let n = perf.len();
-            let dist = outcome_distribution(&perf, &flat(n, eps));
+            let dist = outcome_distribution(&perf, &flat(n, eps)).unwrap();
             let sum: f64 = dist.iter().map(|(_, p)| p).sum();
             assert!(
                 (sum - 1.0).abs() < 1e-6,
@@ -620,7 +649,7 @@ mod tests {
     fn two_team_distribution_matches_the_closed_form() {
         let perf = [g(3.0, 6.0), g(-2.0, 1.0)];
         let eps = 1.5;
-        let dist = outcome_distribution(&perf, &flat(2, eps));
+        let dist = outcome_distribution(&perf, &flat(2, eps)).unwrap();
         let (wa, wb) = closed_form_two(perf[0], perf[1], eps);
 
         let find = |ranks: &[u32]| {
@@ -653,10 +682,10 @@ mod tests {
         let perf = [g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)];
         let eps = 1.5;
         let margins = flat(3, eps);
-        let dist = outcome_distribution(&perf, &margins);
+        let dist = outcome_distribution(&perf, &margins).unwrap();
 
         for (ranks, expected) in &dist {
-            let direct = ranking_probability(&perf, &margins, ranks);
+            let direct = ranking_probability(&perf, &margins, ranks).unwrap();
             assert!(
                 (direct - expected).abs() < 1e-9,
                 "ranks {ranks:?}: direct {direct} vs distribution {expected}"
@@ -675,7 +704,7 @@ mod tests {
         let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)];
         let mut previous = 0.0;
         for eps in [0.0, 0.5, 1.0, 2.0, 4.0, 8.0, 24.0] {
-            let p = ranking_probability(&perf, &flat(3, eps), &[0, 0, 0]);
+            let p = ranking_probability(&perf, &flat(3, eps), &[0, 0, 0]).unwrap();
             assert!(p >= previous, "eps={eps}: {p} < {previous}");
             if eps == 0.0 {
                 assert!(p < 1e-12, "a tie needs a margin, got {p}");
@@ -696,7 +725,7 @@ mod tests {
         let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)];
         let sweep: Vec<f64> = [0.5, 2.0, 4.0, 8.0, 16.0]
             .iter()
-            .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1]))
+            .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1]).unwrap())
             .collect();
         let peak = sweep
             .iter()
@@ -717,7 +746,7 @@ mod tests {
     #[test]
     fn ties_are_impossible_without_a_draw_margin() {
         let perf = [g(0.0, 4.0), g(0.0, 4.0), g(0.0, 4.0)];
-        let dist = outcome_distribution(&perf, &flat(3, 0.0));
+        let dist = outcome_distribution(&perf, &flat(3, 0.0)).unwrap();
         assert_eq!(dist.len(), 6, "expected only the 6 strict orders: {dist:?}");
         assert!(dist.iter().all(|(r, _)| {
             let mut seen = r.clone();

src/rating.rs

@@ -23,7 +23,24 @@ pub struct Rating<T: Time = i64, D: Drift<T> = ConstantDrift> {
 }
 
 impl<T: Time, D: Drift<T>> Rating<T, D> {
+    /// # Panics
+    ///
+    /// Panics unless `beta` is finite and non-negative, matching
+    /// `HistoryBuilder::beta`.
+    ///
+    /// Zero is allowed and meaningful — performance is then exactly skill, and
+    /// the fit differs measurably from a positive beta rather than degenerating.
+    /// Negative is rejected because `beta` enters only as `beta^2`: measured, a
+    /// negative beta returned results **bit identical** to its absolute value,
+    /// and a NaN beta reached `Game::ranked`, which returned `Ok` carrying a
+    /// `Gaussian { pi: NaN, tau: NaN }` — there is no `converge` on that path to
+    /// catch it.
     pub fn new(prior: Gaussian, beta: f64, drift: D) -> Self {
+        assert!(
+            beta.is_finite() && beta >= 0.0,
+            "beta must be finite and non-negative (got {beta}); it is only ever \
+             squared, so a negative value would silently behave as its absolute value"
+        );
         Self {
             prior,
             beta,

tests/cross_process_determinism.rs

@@ -0,0 +1,157 @@
+//! Determinism across *processes*, which an in-process test cannot see.
+//!
+//! Rust seeds its default hasher once per process, so every `HashMap`
+//! iteration order is fixed for a run and varies between runs. A test that
+//! compares results within one process therefore cannot detect a float sum
+//! whose order comes from a map — all its samples share one seed.
+//!
+//! That is not hypothetical. `tests/determinism.rs` compares four thread counts
+//! inside one process and passed throughout, while `posterior_of` was returning
+//! two distinct bit patterns across 40 separate runs on identical input.
+//!
+//! This re-executes the test binary and compares `f64::to_bits`.
+
+use std::{env, process::Command};
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, History, Member, NullObserver, Outcome, Team,
+    UnknownKeys,
+};
+
+/// Set in the child so it reports instead of re-spawning.
+const CHILD: &str = "TSTT_DETERMINISM_CHILD";
+
+const RUNS: usize = 40;
+
+type H = History<i64, ConstantDrift, NullObserver, String>;
+
+fn fitted() -> H {
+    let mut h: H = History::builder_with_key()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.05))
+        .unknown_keys(UnknownKeys::Prior)
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build();
+
+    let mut events = Vec::new();
+    for t in 0..12i64 {
+        for k in 0..6usize {
+            let a = format!("p{}", (t as usize * 6 + k) % 10);
+            let b = format!("p{}", (t as usize * 6 + k + 4) % 10);
+            events.push(Event {
+                time: t,
+                teams: smallvec![
+                    Team::with_members([Member::new(a)]),
+                    Team::with_members([Member::new(b)]),
+                ],
+                outcome: Outcome::scores([3.0, 1.0]),
+            });
+        }
+    }
+    h.add_events(events).unwrap();
+    assert!(h.converge().unwrap().converged);
+    h
+}
+
+/// Every quantity that could plausibly depend on iteration order, as bits.
+fn fingerprint() -> String {
+    let h = fitted();
+
+    // Unknown keys with UNEQUAL but COMPARABLE coefficients, which is what
+    // makes the sum order-sensitive.
+    //
+    // Equal terms sum order-independently and would make this pass vacuously.
+    // Terms of wildly different magnitudes are no better: the small ones fall
+    // below the running total's ULP and are absorbed whatever the order —
+    // measured, spreading these over nine decades dropped the detection rate
+    // to roughly one run in forty. Comparable sizes keep every term able to
+    // change the last bits.
+    let ghosts: Vec<String> = (0..24).map(|i| format!("ghost{i}")).collect();
+    let mut terms: Vec<(&String, f64)> = ghosts
+        .iter()
+        .enumerate()
+        .map(|(i, k)| (k, 1.0 + i as f64 * 0.37))
+        .collect();
+    let known = "p0".to_string();
+    terms.push((&known, -1.0));
+
+    let posterior = h.posterior_of(&terms).unwrap();
+
+    let a = "p0".to_string();
+    let b = "p1".to_string();
+    let target = [(&a, 1.0), (&b, -1.0)];
+    let teams: [&[&String]; 2] = [&[&a], &[&b]];
+    let evr = h.expected_variance_reduction(&teams, &target).unwrap();
+
+    let curves = h.learning_curves();
+    let mut curve_bits: u64 = 0;
+    let mut keys: Vec<&String> = curves.keys().collect();
+    keys.sort();
+    for key in keys {
+        for (t, g) in &curves[key] {
+            curve_bits ^= (*t as u64).rotate_left(17)
+                ^ g.mu().to_bits().rotate_left(31)
+                ^ g.sigma().to_bits();
+        }
+    }
+
+    format!(
+        "post={:016x} evr={:016x} le={:016x} curves={curve_bits:016x}",
+        posterior.sigma().to_bits(),
+        evr.to_bits(),
+        h.log_evidence().to_bits(),
+    )
+}
+
+#[test]
+fn results_are_identical_across_processes() {
+    if env::var(CHILD).is_ok() {
+        println!("FINGERPRINT {}", fingerprint());
+        return;
+    }
+
+    let exe = env::current_exe().expect("current exe");
+    let mut seen: Vec<String> = Vec::new();
+
+    for run in 0..RUNS {
+        let out = Command::new(&exe)
+            .args([
+                "results_are_identical_across_processes",
+                "--exact",
+                "--nocapture",
+            ])
+            .env(CHILD, "1")
+            .output()
+            .expect("spawn child");
+        assert!(
+            out.status.success(),
+            "child {run} failed: {}",
+            String::from_utf8_lossy(&out.stderr)
+        );
+        let stdout = String::from_utf8_lossy(&out.stdout);
+        let line = stdout
+            .lines()
+            .find_map(|l| l.strip_prefix("FINGERPRINT "))
+            .unwrap_or_else(|| panic!("child {run} printed no fingerprint:\n{stdout}"))
+            .to_string();
+        seen.push(line);
+    }
+
+    let first = &seen[0];
+    let differing: Vec<&String> = seen.iter().filter(|s| *s != first).collect();
+    assert!(
+        differing.is_empty(),
+        "results differ across processes on identical input.\n  {} of {RUNS} runs differed\n  \
+         first:     {first}\n  differing: {}",
+        differing.len(),
+        differing[0]
+    );
+}

tests/determinism.rs

@@ -1,101 +1,201 @@
-//! Determinism tests: identical posteriors across RAYON_NUM_THREADS
-//! values. Only compiled with the `rayon` feature.
+//! Determinism across `RAYON_NUM_THREADS`, on a workload that actually reaches
+//! the parallel path.
+//!
+//! This test previously proved less than it appeared to. `sweep_color_groups`
+//! takes its `par_iter` branch only for colour groups of at least
+//! `RAYON_THRESHOLD` (64) events, and the old fixture built 20 slices of 10
+//! events — a colour group is a subset of one slice's events, so it could never
+//! exceed 10. The branch was unreachable, confirmed by CPU-vs-wall time:
+//! `user 0.64` on eight threads is one core.
+//!
+//! It also compared a single competitor's curve out of forty, and never
+//! compared `log_evidence`, `final_step` or `iterations`.
+//!
+//! The fixture below guarantees the parallel branch **by construction**: within
+//! a slice every event uses a disjoint pair of competitors, so greedy colouring
+//! puts all of them in colour 0, and that group is `EVENTS_PER_SLICE` long.
+//! Competitors recur across slices, so the fit still has temporal coupling and
+//! drift rather than being a set of independent duels.
 
 #![cfg(feature = "rayon")]
 
 use smallvec::smallvec;
-use trueskill_tt::{ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team};
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, Gaussian, History, Member, Outcome, Team,
+};
 
-/// Build a deterministic workload using a simple LCG (no external rand crate).
-fn build_and_converge(seed: u64) -> Vec<(i64, trueskill_tt::Gaussian)> {
+/// Comfortably above the crate's internal `RAYON_THRESHOLD` of 64.
+const EVENTS_PER_SLICE: usize = 96;
+const SLICES: i64 = 8;
+/// Two per event, all disjoint within a slice.
+const COMPETITORS: usize = EVENTS_PER_SLICE * 2;
+
+/// Everything a thread count could plausibly perturb.
+struct Fingerprint {
+    curves: Vec<(String, Vec<(i64, Gaussian)>)>,
+    log_evidence: f64,
+    final_step: (f64, f64),
+    iterations: usize,
+}
+
+fn build_and_converge() -> Fingerprint {
     let mut h = History::<i64, _, _, String>::builder_with_key()
         .mu(25.0)
         .sigma(25.0 / 3.0)
         .beta(25.0 / 6.0)
         .drift(ConstantDrift(25.0 / 300.0))
         .convergence(ConvergenceOptions {
-            max_iter: 30,
-            epsilon: 1e-6,
+            max_iter: 20_000,
+            epsilon: 1e-9,
             alpha: 1.0,
         })
         .build();
 
-    // LCG for deterministic pseudo-random ints.
-    let mut rng = seed;
-    let mut next = || {
-        rng = rng
-            .wrapping_mul(6364136223846793005)
-            .wrapping_add(1442695040888963407);
-        rng
-    };
-
-    let mut events: Vec<Event<i64, String>> = Vec::with_capacity(200);
-    for ev_i in 0..200 {
-        let a = (next() % 40) as usize;
-        let mut b = (next() % 40) as usize;
-        while b == a {
-            b = (next() % 40) as usize;
+    let mut events: Vec<Event<i64, String>> = Vec::new();
+    for slice in 0..SLICES {
+        for e in 0..EVENTS_PER_SLICE {
+            // Disjoint within the slice: event `e` owns competitors 2e and
+            // 2e+1. Rotating by the slice index makes the pairings differ
+            // between slices, so competitors accumulate a real history.
+            let a = (2 * e + slice as usize) % COMPETITORS;
+            let b = (2 * e + 1 + slice as usize * 3) % COMPETITORS;
+            if a == b {
+                continue;
+            }
+            events.push(Event {
+                time: slice + 1,
+                teams: smallvec![
+                    Team::with_members([Member::new(format!("p{a}"))]),
+                    Team::with_members([Member::new(format!("p{b}"))]),
+                ],
+                outcome: Outcome::winner(u32::from((e + slice as usize) % 2 == 0), 2),
+            });
         }
-        // ~10 events per slice so color groups have material parallelism.
-        events.push(Event {
-            time: (ev_i as i64 / 10) + 1,
-            teams: smallvec![
-                Team::with_members([Member::new(format!("p{a}"))]),
-                Team::with_members([Member::new(format!("p{b}"))]),
-            ],
-            outcome: Outcome::winner((next() % 2) as u32, 2),
-        });
     }
     h.add_events(events).unwrap();
-    let _ = h.converge().unwrap();
-    // Sample one competitor's curve for the comparison.
-    h.learning_curve("p0")
+
+    let report = h.converge().expect("fixture must converge");
+
+    let mut curves: Vec<(String, Vec<(i64, Gaussian)>)> = h
+        .learning_curves()
+        .into_iter()
+        .map(|(k, v)| (k.clone(), v))
+        .collect();
+    curves.sort_by(|a, b| a.0.cmp(&b.0));
+
+    Fingerprint {
+        curves,
+        log_evidence: h.log_evidence(),
+        final_step: report.final_step,
+        iterations: report.iterations,
+    }
 }
 
 #[test]
 fn posteriors_identical_across_thread_counts() {
     let sizes = [1usize, 2, 4, 8];
-    let mut results: Vec<Vec<(i64, trueskill_tt::Gaussian)>> = Vec::new();
+    let mut results: Vec<Fingerprint> = Vec::new();
+
     for &n in &sizes {
         let pool = rayon::ThreadPoolBuilder::new()
             .num_threads(n)
             .build()
             .expect("rayon pool build");
-        let curve = pool.install(|| build_and_converge(42));
-        results.push(curve);
+        results.push(pool.install(build_and_converge));
     }
 
     let reference = &results[0];
-    for (i, curve) in results.iter().enumerate().skip(1) {
+
+    // Guard against the failure this test previously had: passing while
+    // measuring almost nothing.
+    assert!(
+        reference.curves.len() > 100,
+        "expected every competitor's curve, got {}",
+        reference.curves.len()
+    );
+
+    for (i, got) in results.iter().enumerate().skip(1) {
+        let n = sizes[i];
+
         assert_eq!(
-            curve.len(),
-            reference.len(),
-            "curve length differs at {n} threads",
-            n = sizes[i],
+            got.iterations, reference.iterations,
+            "iterations differ at {n} threads"
         );
-        for (j, (&(t_ref, g_ref), &(t, g))) in reference.iter().zip(curve.iter()).enumerate() {
+        assert_eq!(
+            got.final_step.0.to_bits(),
+            reference.final_step.0.to_bits(),
+            "final_step.0 differs at {n} threads: {:?} vs {:?}",
+            reference.final_step,
+            got.final_step
+        );
+        assert_eq!(
+            got.final_step.1.to_bits(),
+            reference.final_step.1.to_bits(),
+            "final_step.1 differs at {n} threads"
+        );
+        assert_eq!(
+            got.log_evidence.to_bits(),
+            reference.log_evidence.to_bits(),
+            "log_evidence differs at {n} threads: {} vs {}",
+            reference.log_evidence,
+            got.log_evidence
+        );
+
+        assert_eq!(
+            got.curves.len(),
+            reference.curves.len(),
+            "competitor count differs at {n} threads"
+        );
+
+        for ((ref_key, ref_curve), (key, curve)) in reference.curves.iter().zip(got.curves.iter()) {
+            assert_eq!(ref_key, key, "competitor order differs at {n} threads");
             assert_eq!(
-                t_ref,
-                t,
-                "time point {j} differs at {n} threads: ref={t_ref} vs got={t}",
-                n = sizes[i],
-            );
-            assert_eq!(
-                g_ref.mu().to_bits(),
-                g.mu().to_bits(),
-                "mu bits differ at {n} threads, time {t}: ref={ref_mu} got={got_mu}",
-                n = sizes[i],
-                ref_mu = g_ref.mu(),
-                got_mu = g.mu(),
-            );
-            assert_eq!(
-                g_ref.sigma().to_bits(),
-                g.sigma().to_bits(),
-                "sigma bits differ at {n} threads, time {t}: ref={ref_sigma} got={got_sigma}",
-                n = sizes[i],
-                ref_sigma = g_ref.sigma(),
-                got_sigma = g.sigma(),
+                curve.len(),
+                ref_curve.len(),
+                "curve length differs for {key} at {n} threads"
             );
+            for (&(t_ref, g_ref), &(t, g)) in ref_curve.iter().zip(curve.iter()) {
+                assert_eq!(t_ref, t, "time point differs for {key} at {n} threads");
+                assert_eq!(
+                    g_ref.mu().to_bits(),
+                    g.mu().to_bits(),
+                    "mu differs for {key} at t={t}, {n} threads: {} vs {}",
+                    g_ref.mu(),
+                    g.mu()
+                );
+                assert_eq!(
+                    g_ref.sigma().to_bits(),
+                    g.sigma().to_bits(),
+                    "sigma differs for {key} at t={t}, {n} threads: {} vs {}",
+                    g_ref.sigma(),
+                    g.sigma()
+                );
+            }
         }
     }
 }
+
+/// The fixture must keep reaching the parallel branch.
+///
+/// `RAYON_THRESHOLD` is private, so this pins the property that makes the
+/// branch reachable rather than the branch itself: within a slice every event
+/// uses a disjoint competitor pair, so greedy colouring puts all
+/// `EVENTS_PER_SLICE` of them in one colour group. If someone shrinks the
+/// fixture, this fails rather than the suite quietly going back to testing the
+/// sequential path.
+#[test]
+fn the_fixture_still_exceeds_the_rayon_threshold() {
+    const RAYON_THRESHOLD: usize = 64;
+    const {
+        assert!(
+            EVENTS_PER_SLICE >= RAYON_THRESHOLD,
+            "a colour group holds at most EVENTS_PER_SLICE events, which must \
+             reach the crate's RAYON_THRESHOLD for the parallel sweep to run"
+        );
+    }
+
+    // Measured by instrumenting `sweep_color_groups`: this fixture produces
+    // one colour group of 96 events and takes the parallel branch on all 872
+    // sweeps. The old fixture's 10-event slices could not reach 64 at all.
+    assert_eq!(EVENTS_PER_SLICE, 96);
+}

tests/joint_handle.rs

@@ -265,3 +265,75 @@ fn unseen_competitors_match_the_one_shot_path() {
     assert_eq!(one_shot.pi(), cached.pi());
     assert_eq!(one_shot.tau(), cached.tau());
 }
+
+/// A drift too small to represent must collapse, not corrupt the matrix.
+///
+/// The collapse rule used to fire only at `drift <= 0.0` exactly. Anything
+/// smaller-but-positive got an explicit `1.0 / drift` precision, and at
+/// `drift = 1e-16` that entry is `1e16` — so `1e16 + 0.28` rounds back to
+/// `1e16` and the prior and contrasts are annihilated in the stored `f64`.
+///
+/// Measured before the fix, at `drift_scale = 1e-10` this returned a variance
+/// **12 000x too small** (a 111x overconfident interval) as `Ok`, with a band
+/// just above it returning a misleading `JointUnavailable`.
+#[test]
+fn a_drift_too_small_to_represent_collapses_rather_than_corrupting() {
+    fn variance(scale: f64) -> f64 {
+        let mut h: History<i64, ConstantDrift, _, String> = History::builder_with_key()
+            .mu(0.0)
+            .sigma(6.0)
+            .beta(1.0)
+            .score_sigma(2.0)
+            .drift(ConstantDrift(0.5))
+            .convergence(ConvergenceOptions {
+                max_iter: 20_000,
+                epsilon: 1e-13,
+                alpha: 1.0,
+            })
+            .build();
+        let mut events = Vec::new();
+        for t in 0..15i64 {
+            for k in 0..4usize {
+                let x = format!("p{}", (t as usize * 4 + k) % 8);
+                let y = format!("p{}", (t as usize * 4 + k + 3) % 8);
+                events.push(Event {
+                    time: t,
+                    teams: smallvec![
+                        Team::with_members([Member::new(x).with_drift_scale(scale)]),
+                        Team::with_members([Member::new(y).with_drift_scale(scale)]),
+                    ],
+                    outcome: Outcome::scores([3.0, 1.0]),
+                });
+            }
+        }
+        h.add_events(events).unwrap();
+        assert!(h.converge().unwrap().converged);
+        let (a, b) = ("p0".to_string(), "p1".to_string());
+        let joint = h
+            .joint()
+            .expect("a tiny drift must not make the joint unavailable");
+        let g = joint.posterior_of(&[(&a, 1.0), (&b, -1.0)]).unwrap();
+        g.sigma() * g.sigma()
+    }
+
+    let collapsed = variance(0.0);
+
+    // Below the threshold every scale must reach the collapsed answer exactly,
+    // and none may error.
+    for scale in [1e-3, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12] {
+        let v = variance(scale);
+        assert_eq!(
+            v.to_bits(),
+            collapsed.to_bits(),
+            "drift_scale {scale:e}: {v} vs collapsed {collapsed}"
+        );
+    }
+
+    // Above it, real drift is still modelled — otherwise this test would pass
+    // by collapsing everything.
+    let drifting = variance(1e-2);
+    assert!(
+        (drifting - collapsed).abs() / collapsed > 1e-5,
+        "a drift of 1e-2 must still move the answer: {drifting} vs {collapsed}"
+    );
+}

tests/libm_rule.rs

@@ -0,0 +1,175 @@
+//! The libm rule, enforced rather than asserted in prose.
+//!
+//! CLAUDE.md requires transcendentals to go through `libm`, not `std`:
+//!
+//! > IEEE 754 pins the basic operations and `sqrt` but says nothing about
+//! > `exp`/`log`/`erf`, and `std` delegates to the *system* math library —
+//! > measured, `f64::exp` and `libm::exp` disagree on 9.7% of inputs by one
+//! > ULP. Since inference is an iterative fixed point, one ULP can change an
+//! > iteration count.
+//!
+//! The rule was stated clearly and still violated in three production sites,
+//! one of them `hypot` on the path of every scored event — whose measured
+//! divergence, 12.1%, is *higher* than the `exp` figure the rule cites as its
+//! own justification. Prose is evidently not enough, so this is a test.
+//!
+//! Tests may use either, which the crate documents, so `#[cfg(test)]` blocks
+//! are excluded.
+
+use std::{fs, path::Path};
+
+/// Method-call spellings that reach the system math library.
+///
+/// `sqrt` is deliberately absent: IEEE 754 specifies it exactly, so `std` and
+/// `libm` cannot disagree. `abs`, `recip`, `powi` and `mul_add` are likewise
+/// exact or specified.
+const FORBIDDEN: &[&str] = &[
+    "exp", "exp2", "exp_m1", "ln", "ln_1p", "log", "log2", "log10", "powf", "sin", "cos", "tan",
+    "asin", "acos", "atan", "atan2", "sinh", "cosh", "tanh", "hypot", "cbrt", "erf", "erfc",
+];
+
+/// Strip `#[cfg(test)]` items by brace matching, plus comments and string
+/// literals, so a mention in prose is not mistaken for a call.
+fn production_code(source: &str) -> String {
+    let mut out = String::with_capacity(source.len());
+    let bytes: Vec<char> = source.chars().collect();
+    let mut i = 0;
+
+    while i < bytes.len() {
+        let rest: String = bytes[i..].iter().take(16).collect();
+
+        if rest.starts_with("#[cfg(test)]") {
+            // Skip to the opening brace of the guarded item, then past its
+            // matching close.
+            let mut j = i;
+            while j < bytes.len() && bytes[j] != '{' {
+                j += 1;
+            }
+            let mut depth = 0usize;
+            while j < bytes.len() {
+                match bytes[j] {
+                    '{' => depth += 1,
+                    '}' => {
+                        depth -= 1;
+                        if depth == 0 {
+                            j += 1;
+                            break;
+                        }
+                    }
+                    _ => {}
+                }
+                j += 1;
+            }
+            i = j;
+            continue;
+        }
+
+        if rest.starts_with("//") {
+            while i < bytes.len() && bytes[i] != '\n' {
+                i += 1;
+            }
+            continue;
+        }
+
+        if rest.starts_with("/*") {
+            i += 2;
+            while i + 1 < bytes.len() && !(bytes[i] == '*' && bytes[i + 1] == '/') {
+                i += 1;
+            }
+            i += 2;
+            continue;
+        }
+
+        if bytes[i] == '"' {
+            i += 1;
+            while i < bytes.len() && bytes[i] != '"' {
+                if bytes[i] == '\\' {
+                    i += 1;
+                }
+                i += 1;
+            }
+            i += 1;
+            continue;
+        }
+
+        out.push(bytes[i]);
+        i += 1;
+    }
+
+    out
+}
+
+fn rust_files(dir: &Path, out: &mut Vec<std::path::PathBuf>) {
+    for entry in fs::read_dir(dir).expect("read src") {
+        let path = entry.expect("dir entry").path();
+        if path.is_dir() {
+            rust_files(&path, out);
+        } else if path.extension().is_some_and(|e| e == "rs") {
+            out.push(path);
+        }
+    }
+}
+
+#[test]
+fn production_code_never_calls_a_std_transcendental() {
+    let mut files = Vec::new();
+    rust_files(Path::new("src"), &mut files);
+    assert!(files.len() > 10, "expected to find the crate's sources");
+
+    let mut offences = Vec::new();
+
+    for path in &files {
+        let source = fs::read_to_string(path).expect("read source");
+        let code = production_code(&source);
+
+        for (n, line) in code.lines().enumerate() {
+            for name in FORBIDDEN {
+                let needle = format!(".{name}(");
+                if line.contains(&needle) {
+                    offences.push(format!("{}:{}: {}", path.display(), n + 1, line.trim()));
+                }
+            }
+        }
+    }
+
+    assert!(
+        offences.is_empty(),
+        "production code must call libm, not std, for transcendentals \
+         (`sqrt` is exempt — IEEE 754 specifies it):\n{}",
+        offences.join("\n")
+    );
+}
+
+/// The stripper has to actually strip, or the test above passes vacuously.
+#[test]
+fn the_test_module_stripper_works() {
+    let source = r#"
+fn production() { let _ = libm::exp(1.0); }
+
+#[cfg(test)]
+mod tests {
+    fn allowed() { let x = 1.0f64.exp(); }
+}
+
+fn also_production() {}
+"#;
+    let code = production_code(source);
+    assert!(
+        code.contains("also_production"),
+        "stripped too much: {code}"
+    );
+    assert!(
+        !code.contains(".exp()"),
+        "failed to strip cfg(test): {code}"
+    );
+}
+
+/// And it must not strip a doc comment's worth of prose into oblivion, nor
+/// mistake prose for a call.
+#[test]
+fn prose_is_not_mistaken_for_a_call() {
+    let source = "/// Uses `x.exp()` in the docs.\nfn f() { let _ = libm::exp(1.0); }\n";
+    let code = production_code(source);
+    assert!(!code.contains(".exp()"), "doc comment leaked: {code}");
+    assert!(code.contains("libm::exp"), "stripped real code: {code}");
+}

tests/non_finite_results.rs

@@ -10,7 +10,9 @@
 //! `!tuple_gt(..)`. These tests pin the guard from outside.
 
 use smallvec::smallvec;
-use trueskill_tt::{Event, Gaussian, History, InferenceError, Member, Outcome, Team};
+use trueskill_tt::{
+    ConstantDrift, Event, Gaussian, History, InferenceError, Member, Outcome, Team,
+};
 
 fn scored_fit(
     sigma: f64,
@@ -115,3 +117,101 @@ fn merely_extreme_parameters_still_converge() {
     assert!(scored_fit(6.0, 1.0, 1e6, [3.0, 1.0]).unwrap());
     assert!(scored_fit(6.0, 1.0, 1.0, [1e150, -1e150]).unwrap());
 }
+
+/// A NaN in one competitor must not be masked by a healthy competitor reduced
+/// after it.
+///
+/// The convergence step is a fold over a `HashMap`, so which competitor is
+/// reduced last is per-process hash order. Before the fix, `tuple_max` dropped
+/// a NaN accumulator in favour of the next finite delta and this returned
+/// `Ok(converged: true)` with a NaN posterior in **16 of 30 runs** on identical
+/// input. Deterministic now, but note this test can only ever sample one hash
+/// order per run — the ordering guarantee itself is pinned by
+/// `tuple_max_propagates_a_nan_from_any_position` in the crate's unit tests.
+#[test]
+fn a_nan_competitor_is_not_masked_by_a_healthy_one() {
+    let mut h = History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .p_draw(0.1)
+        .build();
+    h.add_events(vec![
+        Event {
+            time: 1i64,
+            teams: smallvec![
+                Team::with_members([Member::new("a").with_prior(Gaussian::from_ms(0.0, 1e-200))]),
+                Team::with_members([Member::new("b")]),
+            ],
+            outcome: Outcome::winner(0, 2),
+        },
+        // A healthy pair in the same slice, to be reduced alongside the NaN.
+        Event {
+            time: 1i64,
+            teams: smallvec![
+                Team::with_members([Member::new("c")]),
+                Team::with_members([Member::new("d")]),
+            ],
+            outcome: Outcome::winner(0, 2),
+        },
+    ])
+    .unwrap();
+
+    let err = h
+        .converge()
+        .expect_err("a NaN fit must never be reported as converged");
+    assert!(
+        matches!(err, InferenceError::NonFiniteResult { .. }),
+        "{err:?}"
+    );
+}
+
+/// A tie observed with a narrow draw margin between far-apart competitors must
+/// produce a fit, not NaN skills.
+///
+/// The tie branch forms the truncated variance from `v^2 - u`, and both grow as
+/// `alpha^2` while their difference stays `O(1)`. Deep enough into the tail
+/// that subtraction had four digits left: measured, it returned `1 - w`
+/// negative and `sqrt` of it was NaN. The half-line escape hatch did not cover
+/// it, because that keys on how many window-widths from the mean the window
+/// sits and a narrow window fails that however deep it is.
+///
+/// These parameters are ordinary for a precise-scoring domain, and the
+/// neighbouring wider-margin case always worked — so this was a cliff, not
+/// "extreme inputs break".
+#[test]
+fn a_narrow_draw_margin_far_into_the_tail_still_fits() {
+    for (beta, p_draw, sd, gap) in [
+        (1e-2, 1e-8, 1e-2, 10.0),
+        (1e-3, 1e-9, 1e-3, 1.0),
+        (1e-4, 1e-12, 1e-4, 1.0),
+    ] {
+        let mut h = History::builder()
+            .mu(0.0)
+            .sigma(sd)
+            .beta(beta)
+            .p_draw(p_draw)
+            .drift(ConstantDrift(0.0))
+            .build();
+        h.add_events(vec![Event {
+            time: 1i64,
+            teams: smallvec![
+                Team::with_members([Member::new("a").with_prior(Gaussian::from_ms(0.0, sd))]),
+                Team::with_members([Member::new("b").with_prior(Gaussian::from_ms(gap, sd))]),
+            ],
+            outcome: Outcome::draw(2),
+        }])
+        .unwrap();
+
+        let report = h
+            .converge()
+            .unwrap_or_else(|e| panic!("beta {beta:e}, p_draw {p_draw:e}: {e:?}"));
+        assert!(report.converged);
+
+        let skill = h.current_skill(&"a").unwrap();
+        assert!(
+            skill.mu().is_finite() && skill.sigma().is_finite() && skill.sigma() > 0.0,
+            "beta {beta:e}, p_draw {p_draw:e}: {skill:?}"
+        );
+    }
+}

tests/prediction_bounds.rs

@@ -0,0 +1,154 @@
+//! Bounds that any correct implementation must satisfy, swept rather than
+//! spot-checked.
+//!
+//! The crate's docs call the `ln k` ceiling "the sharpest available test of an
+//! implementation", and record that an early prototype returned 4.77 nats. It
+//! was violated again — 3.237828 nats against `ln 2` — because the existing
+//! check sampled one fixture and the violation lives in a specific regime: a
+//! large ratio between the widest and narrowest performance sigma, where the
+//! shared prediction grid could not resolve the narrow density and returned
+//! probabilities greater than one.
+//!
+//! A single fixture cannot defend a bound like this. A sweep can.
+
+use trueskill_tt::{
+    ConstantDrift, GameOptions, Gaussian, InferenceError, Rating, expected_information_gain,
+};
+
+type R = Rating<i64, ConstantDrift>;
+
+/// How many random matchups the ceiling sweep draws.
+///
+/// Scaled by build profile rather than fixed. Each sample runs a full inference
+/// pass per outcome, and that is about **19x** faster in release — measured,
+/// 20 000 samples take 12.1s released against 23s for 2 000 in debug. `just
+/// test` runs three debug feature combinations and one release one, so a fixed
+/// count pays the slow price three times and the fast one once, which is
+/// exactly backwards.
+///
+/// The debug run is here to prove the sweep still compiles and holds on a small
+/// sample; the release run is the one that actually searches. The violation
+/// this guards was found at a rate near 1.8%, so even the debug count expects
+/// tens of hits in the regime.
+#[cfg(debug_assertions)]
+const SAMPLES: usize = 1_000;
+#[cfg(not(debug_assertions))]
+const SAMPLES: usize = 50_000;
+
+/// Deterministic LCG, so a failure is reproducible from the printed seed.
+struct Lcg(u64);
+
+impl Lcg {
+    fn next_f64(&mut self) -> f64 {
+        self.0 = self
+            .0
+            .wrapping_mul(6_364_136_223_846_793_005)
+            .wrapping_add(1_442_695_040_888_963_407);
+        // Top 53 bits to [0, 1).
+        ((self.0 >> 11) as f64) / ((1u64 << 53) as f64)
+    }
+
+    fn in_range(&mut self, lo: f64, hi: f64) -> f64 {
+        lo + (hi - lo) * self.next_f64()
+    }
+
+    /// Log-uniform, so the sweep spends its samples across magnitudes rather
+    /// than crowding the top of the range — the violations live at small sigma.
+    fn log_uniform(&mut self, lo: f64, hi: f64) -> f64 {
+        let t = self.next_f64();
+        (lo.ln() + t * (hi.ln() - lo.ln())).exp()
+    }
+}
+
+#[test]
+fn information_gain_never_exceeds_the_entropy_of_the_outcome() {
+    let mut rng = Lcg(0x5eed_1234_abcd_ef01);
+    let ceiling = 2.0_f64.ln();
+    let mut evaluated = 0usize;
+    let mut refused = 0usize;
+
+    for i in 0..SAMPLES {
+        let mu_a = rng.in_range(-100.0, 100.0);
+        let mu_b = rng.in_range(-100.0, 100.0);
+        let sigma_a = rng.log_uniform(1e-4, 1e2);
+        let sigma_b = rng.log_uniform(1e-4, 1e2);
+        let beta = rng.log_uniform(1e-4, 1e1);
+
+        let a = R::new(Gaussian::from_ms(mu_a, sigma_a), beta, ConstantDrift(0.0));
+        let b = R::new(Gaussian::from_ms(mu_b, sigma_b), beta, ConstantDrift(0.0));
+        let options = GameOptions {
+            p_draw: 0.0,
+            ..GameOptions::default()
+        };
+
+        match expected_information_gain(&[&[a], &[b]], &options) {
+            Ok(gain) => {
+                evaluated += 1;
+                assert!(
+                    gain.is_finite(),
+                    "sample {i}: non-finite gain {gain} \
+                     (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})"
+                );
+                assert!(
+                    gain >= 0.0,
+                    "sample {i}: negative gain {gain} \
+                     (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})"
+                );
+                assert!(
+                    gain <= ceiling + 1e-9,
+                    "sample {i}: gain {gain} exceeds ln 2 = {ceiling} \
+                     (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})"
+                );
+            }
+            // Refusing to answer is acceptable; answering wrongly is not.
+            Err(InferenceError::GridTooCoarse { .. }) => refused += 1,
+            Err(e) => panic!("sample {i}: unexpected error {e:?}"),
+        }
+    }
+
+    // The sweep must actually exercise the function, not pass by refusing
+    // everything.
+    assert!(
+        evaluated * 2 > SAMPLES,
+        "only {evaluated} of {SAMPLES} samples were evaluated ({refused} refused); \
+         the sweep is no longer testing anything"
+    );
+    // And it must still reach the regime where the ceiling was violated —
+    // large sigma ratios, which is exactly where the grid now refuses. Without
+    // this the sweep could drift into only-easy inputs and stop being a guard.
+    assert!(
+        refused > 0,
+        "no sample reached the coarse-grid regime; the sweep no longer covers \
+         the case that produced 3.24 nats"
+    );
+}
+
+/// The regime that produced 3.237828 nats, pinned exactly.
+#[test]
+fn the_known_ceiling_violation_no_longer_answers_wrongly() {
+    let a = R::new(
+        Gaussian::from_ms(9.577_887_112_129_012, 0.000_132_507_526_585_134_38),
+        0.000_307_235_559_013_096_2,
+        ConstantDrift(0.0),
+    );
+    let b = R::new(
+        Gaussian::from_ms(-14.114_932_828_525_696, 91.586_690_140_921_16),
+        0.000_307_235_559_013_096_2,
+        ConstantDrift(0.0),
+    );
+    let options = GameOptions {
+        p_draw: 0.0,
+        ..GameOptions::default()
+    };
+
+    match expected_information_gain(&[&[a], &[b]], &options) {
+        Ok(gain) => assert!(
+            gain <= 2.0_f64.ln() + 1e-9,
+            "returned {gain}, over the ln 2 ceiling"
+        ),
+        Err(InferenceError::GridTooCoarse { needed, max }) => {
+            assert!(needed > max, "needed {needed} should exceed max {max}");
+        }
+        Err(e) => panic!("unexpected error {e:?}"),
+    }
+}

tests/quality.rs

@@ -164,3 +164,52 @@ fn quality_matches_the_reference_implementation() {
     let refs: Vec<&[Gaussian]> = five.iter().map(Vec::as_slice).collect();
     assert!((quality(&refs, beta) - 0.040).abs() < 1e-9);
 }
+
+/// `quality()` used to compute `det(ata) / det(middle)` in linear space. Both
+/// are products of `k - 1` diagonal entries, so they leave `f64`'s range long
+/// before their ratio does — and the ratio is the only thing the answer needs.
+///
+/// Measured before the fix: at the crate defaults 150 groups was correct, 200
+/// returned `0`, and 250 returned `NaN` where the truth is `9.51e-88`. With a
+/// small beta it bit sooner — `sigma = beta = 1e-3` returned `NaN` at 60 groups
+/// against a true `1.32e-9`, a value that is entirely ordinary.
+///
+/// For `k` single-member groups with equal means the answer has a closed form,
+/// `(beta / sqrt(beta^2 + sigma^2))^(k-1)`, so this checks against arithmetic
+/// rather than against a recorded output.
+#[test]
+fn quality_matches_its_closed_form_past_the_overflow_point() {
+    for (sigma, beta) in [(25.0 / 3.0, 25.0 / 6.0), (1e-3, 1e-3), (50.0, 25.0 / 6.0)] {
+        let rating = vec![Gaussian::from_ms(25.0, sigma)];
+        for k in [2usize, 50, 60, 150, 200, 250, 300] {
+            let groups: Vec<&[Gaussian]> = (0..k).map(|_| rating.as_slice()).collect();
+            let got = quality(&groups, beta);
+            let expected = (beta / (beta * beta + sigma * sigma).sqrt()).powi(k as i32 - 1);
+
+            assert!(
+                got.is_finite(),
+                "sigma {sigma}, beta {beta}, {k} groups: got {got}"
+            );
+            // Subnormal results have no relative precision left to check.
+            if expected > f64::MIN_POSITIVE {
+                let rel = ((got - expected) / expected).abs();
+                assert!(
+                    rel < 1e-11,
+                    "sigma {sigma}, beta {beta}, {k} groups: got {got:e}, \
+                     closed form {expected:e}, rel {rel:e}"
+                );
+            }
+        }
+    }
+}
+
+/// The overflow was in the intermediates, never in the answer: every value
+/// above is an ordinary float. This pins the specific case that returned `NaN`
+/// where the true answer is nine orders of magnitude inside the normal range.
+#[test]
+fn a_small_beta_does_not_overflow_at_sixty_groups() {
+    let rating = vec![Gaussian::from_ms(25.0, 1e-3)];
+    let groups: Vec<&[Gaussian]> = (0..60).map(|_| rating.as_slice()).collect();
+    let got = quality(&groups, 1e-3);
+    assert!((got - 1.317_089e-9).abs() / 1.317_089e-9 < 1e-6, "{got:e}");
+}

tests/validation.rs

@@ -255,3 +255,88 @@ mod builder_parameters {
         );
     }
 }
+
+/// The constructors below `HistoryBuilder`, which 0.8.0's validation did not
+/// reach.
+///
+/// `sigma`, `beta` and `gamma` all enter inference only as squares, so a
+/// negative value behaves as its absolute value and the sign vanishes without
+/// comment. Measured before these guards: `from_ms(25.0, -8.33)` and
+/// `Rating::new(_, -4.17, _)` returned results bit identical to their positive
+/// counterparts, and `Rating::new(_, NaN, _)` reached `Game::ranked`, which
+/// returned `Ok` carrying `Gaussian { pi: NaN, tau: NaN }`.
+mod constructor_parameters {
+    use trueskill_tt::{ConstantDrift, Gaussian, History, InferenceError, Rating};
+
+    #[test]
+    #[should_panic(expected = "sigma must not be negative")]
+    fn a_negative_sigma_is_rejected_by_from_ms() {
+        let _ = Gaussian::from_ms(25.0, -8.33);
+    }
+
+    /// NaN must pass, and that is deliberate: a broken fit produces a NaN
+    /// sigma and `converge` reports it as `NonFiniteResult`. Rejecting it here
+    /// would turn reporting into a panic inside inference.
+    #[test]
+    fn a_nan_sigma_passes_through_from_ms() {
+        let g = Gaussian::from_ms(25.0, f64::NAN);
+        assert!(g.sigma().is_nan() || g.pi().is_nan());
+    }
+
+    #[test]
+    #[should_panic(expected = "beta must be finite and non-negative")]
+    fn a_negative_beta_is_rejected_by_rating_new() {
+        let _ = Rating::<i64, ConstantDrift>::new(Gaussian::default(), -4.17, ConstantDrift(0.0));
+    }
+
+    #[test]
+    #[should_panic(expected = "beta must be finite and non-negative")]
+    fn a_nan_beta_is_rejected_by_rating_new() {
+        let _ =
+            Rating::<i64, ConstantDrift>::new(Gaussian::default(), f64::NAN, ConstantDrift(0.0));
+    }
+
+    #[test]
+    fn a_zero_beta_is_accepted_by_rating_new() {
+        let _ = Rating::<i64, ConstantDrift>::new(Gaussian::default(), 0.0, ConstantDrift(0.0));
+    }
+
+    /// `HistoryBuilder::drift` is generic and cannot inspect an arbitrary
+    /// `Drift`, so the check is on the variance each competitor actually
+    /// accumulates. That also covers a custom implementation.
+    #[test]
+    fn a_non_finite_drift_is_rejected_at_convergence() {
+        for gamma in [f64::NAN, f64::INFINITY] {
+            let mut h = History::builder()
+                .mu(25.0)
+                .sigma(25.0 / 3.0)
+                .beta(25.0 / 6.0)
+                .drift(ConstantDrift(gamma))
+                .build();
+            h.record_winner(&"a", &"b", 1).unwrap();
+            h.record_winner(&"a", &"b", 5).unwrap();
+            let err = h.converge().unwrap_err();
+            assert!(
+                matches!(
+                    err,
+                    InferenceError::InvalidParameter {
+                        name: "drift variance",
+                        ..
+                    }
+                ),
+                "gamma {gamma}: {err:?}"
+            );
+        }
+    }
+
+    /// An ordinary drift is untouched.
+    #[test]
+    fn an_ordinary_drift_still_converges() {
+        let mut h = History::builder()
+            .drift(ConstantDrift(25.0 / 300.0))
+            .build();
+        h.record_winner(&"a", &"b", 1).unwrap();
+        h.record_winner(&"a", &"b", 5).unwrap();
+        assert!(h.converge().unwrap().converged);
+    }
+}