chore: Release trueskill-tt version 0.4.2

Scouting crates.io for a more accurate `erfc` than `libm` turned up a
1.16e-12 relative error in `erfc_inv` at `p_draw = 0.999999`, measured
against a 70-digit `decimal` reference. It looked like too few Newton
steps. It is not: adding a fourth changed nothing.

The error is in forming the argument. `1.0 - 0.999999` is
`1.0000000000287557e-06` — 0.999999 is not representable, and
subtracting from one cancels, leaving 2.9e-11 of relative error before
`erfc_inv` is entered. Given an exactly-representable argument it
returns 1.8e-16. So the routine was never the problem, and the extra
iteration has been reverted rather than shipped as a fix for a defect
that was not there.

`puruspe::inverfc` returns the identical wrong value for the identical
reason, which is what makes the shared upstream cause obvious.

Adds a test that separates the two, and corrects a quantile constant in
`erfc_inv_matches_known_quantiles` that was recalled rather than
computed: `Phi^-1(0.9999995)` is 4.89163847569859, not
4.891638475699099. The others were checked against the same reference
and were right.

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

CHANGELOG.md

@@ -2,6 +2,31 @@
 
 All notable changes to this project will be documented in this file.
 
+## 0.4.2 - 2026-09-07
+
+### Bug Fixes
+
+- fix: replace the erfc approximation with libm, for free
+- fix: route every transcendental through libm, and combine sigmas with hypot
+
+### Testing
+
+- test: localise the erfc_inv tail residual to the caller's argument
+
+## 0.4.1 - 2026-09-07
+
+### Bug Fixes
+
+- fix: correct erfc_inv's sign error and keep evidence in log space
+
+### Miscellaneous Tasks
+
+- chore: Release trueskill-tt version 0.4.1
+
+### Testing
+
+- test: pin quality()'s N-group closed form, closing the README cross-check
+
 ## 0.4.0 - 2026-09-07
 
 ### Breaking Changes
@@ -25,6 +50,10 @@ All notable changes to this project will be documented in this file.
 - feat: add expected information gain for active matchup selection
 - feat: let observers be shared, boxed, or borrowed
 
+### Miscellaneous Tasks
+
+- chore: Release trueskill-tt version 0.4.0
+
 ## 0.3.0 - 2026-09-01
 
 ### Breaking Changes

CLAUDE.md

@@ -97,6 +97,12 @@ History  →  TimeSlice[]  →  Event[]  →  Item[]
   chain underflows to zero, and `ln(0)` is `-inf`.
 - **Colors are contiguous.** `recompute_color_groups` reorders events so each
   color occupies one range; `ColorGroups::groups_are_contiguous` asserts it.
+- **Transcendentals 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. Use `libm::exp` / `libm::log` in
+  inference code; `f64::sqrt` is fine (IEEE specifies it). Tests may use either.
 - **The crate is `#![forbid(unsafe_code)]`.** Keep it that way.
 - **Ingestion order must not change the answer.** Events added one at a time
   must converge to the same fixed point as the same events batched — see

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "trueskill-tt"
-version = "0.4.0"
+version = "0.4.2"
 edition = "2024"
 rust-version = "1.85"
 description = "TrueSkill Through Time: Bayesian skill rating that tracks how skill evolves over time, via Gaussian message passing"
@@ -51,6 +51,7 @@ harness = false
 
 [dependencies]
 approx = { version = "0.5.1", optional = true }
+libm = "0.2.16"
 rayon = { version = "1", optional = true }
 smallvec = "1"
 

README.md

@@ -242,7 +242,7 @@ expensive than `quality()`. Scoring every pairing among `n` competitors is
 - [x] Add Observer (`Observer` / `NullObserver`)
 - [x] Benchmark the inference loop (`benches/batch.rs`, `benches/history_converge.rs`, `benches/ingest.rs`)
 - [x] N-team `predict_outcome` with draw mass, and `expected_information_gain`
-- [ ] Cross-check `quality()` against [sublee/trueskill](https://github.com/sublee/trueskill/tree/master) — N-group support works and is covered by invariants, but no reference values are asserted
+- [x] Cross-check `quality()` against [sublee/trueskill](https://github.com/sublee/trueskill/tree/master) — N identical teams follow the closed form `(1/5)^((n-1)/2)` for the conventional parameters, asserted for n = 2..10, and the n=3/n=5 values (0.200, 0.040) match the reference package
 
 ## License
 

src/acquisition.rs

@@ -47,7 +47,7 @@ fn kl_divergence(q: Gaussian, p: Gaussian) -> f64 {
     }
 
     let mean_gap = q.mu() - p.mu();
-    0.5 * ((var_p / var_q).ln() + (var_q + mean_gap * mean_gap) / var_p - 1.0)
+    0.5 * (libm::log(var_p / var_q) + (var_q + mean_gap * mean_gap) / var_p - 1.0)
 }
 
 /// Expected information gain of a hypothetical matchup, in nats.

src/factor/margin.rs

@@ -2,7 +2,7 @@ use crate::{
     N_INF,
     factor::{Factor, VarId, VarStore},
     gaussian::Gaussian,
-    pdf,
+    ln_pdf,
 };
 
 /// Gaussian observation factor on a diff variable.
@@ -16,7 +16,7 @@ pub struct MarginFactor {
     pub m_obs: f64,
     pub sigma: f64,
     pub(crate) msg: Gaussian,
-    pub(crate) evidence_cached: Option<f64>,
+    pub(crate) log_evidence_cached: Option<f64>,
 }
 
 impl MarginFactor {
@@ -28,7 +28,7 @@ impl MarginFactor {
             m_obs,
             sigma,
             msg: N_INF,
-            evidence_cached: None,
+            log_evidence_cached: None,
         }
     }
 }
@@ -41,8 +41,8 @@ impl MarginFactor {
         let marginal = vars.get(self.diff);
         let cavity = marginal / self.msg;
 
-        if self.evidence_cached.is_none() {
-            self.evidence_cached = Some(cavity_evidence(cavity, self.m_obs, self.sigma));
+        if self.log_evidence_cached.is_none() {
+            self.log_evidence_cached = Some(cavity_log_evidence(cavity, self.m_obs, self.sigma));
         }
 
         let new_msg = Gaussian::from_ms(self.m_obs, self.sigma);
@@ -61,17 +61,32 @@ impl Factor for MarginFactor {
     }
 
     fn log_evidence(&self, _vars: &VarStore) -> f64 {
-        self.evidence_cached.unwrap_or(1.0).ln()
+        self.log_evidence_cached.unwrap_or(0.0)
     }
 }
 
-/// Density of the observed margin under the cavity, clamped to a positive
-/// floor so a far-out observation cannot underflow to `0.0` and make
-/// `log_evidence` `-inf`.
-fn cavity_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 {
-    let combined_sigma = (cavity.sigma().powi(2) + sigma.powi(2)).sqrt();
+/// `ln` of the observed margin's density under the cavity.
+///
+/// Computed in log space rather than as `pdf(..).ln()`. The density underflows
+/// to zero past about 38 sigma of separation, and clamping that to
+/// `f64::MIN_POSITIVE` reported -708 nats however far out the observation
+/// actually was — 4292 nats adrift at 100 sigma, and unbounded beyond. A score
+/// far from what the model expected is exactly the observation a log-evidence
+/// figure exists to notice.
+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 value = ln_pdf(m_obs, cavity.mu(), combined_sigma);
 
-    pdf(m_obs, cavity.mu(), combined_sigma).max(f64::MIN_POSITIVE)
+    // A degenerate cavity (infinite sigma) is the only way to reach a
+    // non-finite result; fall back to the old floor rather than emit -inf.
+    if value.is_finite() {
+        value
+    } else {
+        f64::MIN_POSITIVE.ln()
+    }
 }
 
 #[cfg(test)]
@@ -113,16 +128,16 @@ mod tests {
         let mut vars = VarStore::new();
         let diff = vars.alloc(Gaussian::from_ms(0.0, 6.0));
         let mut f = MarginFactor::new(diff, 5.0, 1.0);
-        assert!(f.evidence_cached.is_none());
+        assert!(f.log_evidence_cached.is_none());
 
         f.propagate(&mut vars);
-        let z = f.evidence_cached.unwrap();
-        // pdf(5, 0, sqrt(37)) ≈ 0.046783
-        assert!((z - 0.04678300292616668).abs() < 1e-10);
+        let z = f.log_evidence_cached.unwrap();
+        // ln pdf(5, 0, sqrt(37)) = ln(0.046783...)
+        assert!((z.exp() - 0.04678300292616668).abs() < 1e-10);
 
         // Subsequent propagations don't change it.
         f.propagate(&mut vars);
-        assert_eq!(f.evidence_cached.unwrap(), z);
+        assert_eq!(f.log_evidence_cached.unwrap(), z);
     }
 
     #[test]

src/factor/trunc.rs

@@ -1,8 +1,8 @@
 use crate::{
-    N_INF, approx, cdf,
+    N_INF, approx,
     factor::{Factor, VarId, VarStore},
     gaussian::Gaussian,
-    sf,
+    ln_interval, ln_sf,
 };
 
 /// EP truncation factor on a diff variable.
@@ -19,7 +19,7 @@ pub struct TruncFactor {
     /// Outgoing message to the diff variable (initial: `N_INF`, the EP identity).
     pub(crate) msg: Gaussian,
     /// Cached evidence (linear, not log) computed from the cavity on first propagation.
-    pub(crate) evidence_cached: Option<f64>,
+    pub(crate) log_evidence_cached: Option<f64>,
 }
 
 impl TruncFactor {
@@ -30,7 +30,7 @@ impl TruncFactor {
             margin,
             tie,
             msg: N_INF,
-            evidence_cached: None,
+            log_evidence_cached: None,
         }
     }
 }
@@ -43,8 +43,8 @@ impl TruncFactor {
         let marginal = vars.get(self.diff);
         let cavity = marginal / self.msg;
 
-        if self.evidence_cached.is_none() {
-            self.evidence_cached = Some(cavity_evidence(cavity, self.margin, self.tie));
+        if self.log_evidence_cached.is_none() {
+            self.log_evidence_cached = Some(cavity_log_evidence(cavity, self.margin, self.tie));
         }
 
         let trunc = approx(cavity, self.margin, self.tie);
@@ -69,38 +69,34 @@ impl Factor for TruncFactor {
     }
 
     fn log_evidence(&self, _vars: &VarStore) -> f64 {
-        self.evidence_cached.unwrap_or(1.0).ln()
+        self.log_evidence_cached.unwrap_or(0.0)
     }
 }
 
-/// P(diff > margin) for non-tie, P(|diff| < margin) for tie.
+/// `ln P(diff > margin)` for a win, `ln P(|diff| < margin)` for a tie.
 ///
-/// Both branches pick whichever tail keeps their terms *small*, because the
-/// alternative is subtracting two numbers that both approach 1. That
-/// subtraction is not a rounding detail: it loses every digit of an unlikely
-/// outcome's evidence, and an unlikely outcome is precisely the one worth
-/// scoring. `1 - cdf` returned exactly zero past ~8.3 sigma, where the true
-/// probability is 1e-19; clamped, that reached `log_evidence` as -708 instead
-/// of -43.
-///
-/// The clamp remains as a guard rather than a workaround: `erfc` carries ~1e-7
-/// relative error, so a probability of exactly 1 can still come back a hair
-/// above it, and `ln` of a negative would poison the sum for the whole history.
-fn cavity_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 {
+/// Computed in log space throughout. Two earlier shapes both lost the tail:
+/// `1 - cdf(..)` cancelled away every digit of an unlikely outcome, and even
+/// once that was fixed the linear probability underflows to zero past about 38
+/// sigma, where clamping reported -708 nats regardless of the truth. An upset
+/// is the observation a log-evidence figure exists to notice, so it has to stay
+/// exact precisely where it is smallest.
+fn cavity_log_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 {
     let (mu, sigma) = (diff.mu(), diff.sigma());
 
-    let raw = if tie {
-        if mu < -margin {
-            // Both CDFs sit against 1 here; both survival terms are small.
-            sf(-margin, mu, sigma) - sf(margin, mu, sigma)
+    let value = if tie {
+        ln_interval(-margin, margin, mu, sigma)
     } else {
-            cdf(margin, mu, sigma) - cdf(-margin, mu, sigma)
-        }
-    } else {
-        sf(margin, mu, sigma)
+        ln_sf(margin, mu, sigma)
     };
 
-    raw.clamp(f64::MIN_POSITIVE, 1.0)
+    // A degenerate cavity is the only route to a non-finite result; keep the
+    // old floor for it rather than letting -inf poison the whole history's sum.
+    if value.is_finite() {
+        value
+    } else {
+        f64::MIN_POSITIVE.ln()
+    }
 }
 
 #[cfg(test)]
@@ -131,19 +127,19 @@ mod tests {
         let diff = vars.alloc(Gaussian::from_ms(2.0, 3.0));
 
         let mut f = TruncFactor::new(diff, 0.0, false);
-        assert!(f.evidence_cached.is_none());
+        assert!(f.log_evidence_cached.is_none());
 
         f.propagate(&mut vars);
-        assert!(f.evidence_cached.is_some());
-        let first = f.evidence_cached.unwrap();
+        assert!(f.log_evidence_cached.is_some());
+        let first = f.log_evidence_cached.unwrap();
 
         // Evidence should be P(diff > 0) for diff ~ N(2, 9) ≈ 0.748
-        assert!(first > 0.7);
-        assert!(first < 0.8);
+        assert!(first.exp() > 0.7);
+        assert!(first.exp() < 0.8);
 
         // Subsequent propagations don't change it.
         f.propagate(&mut vars);
-        assert_eq!(f.evidence_cached.unwrap(), first);
+        assert_eq!(f.log_evidence_cached.unwrap(), first);
     }
 
     /// The defect this guards: `1 - cdf` collapsed to zero for a surprising
@@ -154,7 +150,7 @@ mod tests {
     #[test]
     fn evidence_of_an_upset_is_not_flattened_to_the_clamp_floor() {
         // diff ~ N(-9, 1) with margin 0: the favoured side lost by nine sigma.
-        let evidence = cavity_evidence(Gaussian::from_ms(-9.0, 1.0), 0.0, false);
+        let evidence = cavity_log_evidence(Gaussian::from_ms(-9.0, 1.0), 0.0, false).exp();
 
         assert!(
             evidence > f64::MIN_POSITIVE,
@@ -175,23 +171,48 @@ mod tests {
     /// Evidence must stay finite and positive however extreme the mismatch,
     /// since `log_evidence` sums across the whole history and one `-inf` or
     /// `NaN` poisons all of it.
+    ///
+    /// Finiteness alone is too weak a bar — the clamped version was finite too,
+    /// and wrong by hundreds of nats. `log_evidence_tracks_the_analytic_tail`
+    /// below is the assertion that actually holds this up.
     #[test]
     fn evidence_stays_positive_and_finite_at_any_separation() {
         for mu in [-300.0f64, -50.0, -9.0, 0.0, 9.0, 50.0, 300.0] {
             for tie in [false, true] {
-                let e = cavity_evidence(Gaussian::from_ms(mu, 1.0), 1.0, tie);
+                let ln_e = cavity_log_evidence(Gaussian::from_ms(mu, 1.0), 1.0, tie);
                 assert!(
-                    e.is_finite() && e > 0.0 && e <= 1.0,
-                    "mu={mu} tie={tie}: evidence {e} is not a probability"
-                );
-                assert!(
-                    e.ln().is_finite(),
-                    "mu={mu} tie={tie}: ln evidence is not finite"
+                    ln_e.is_finite() && ln_e <= 0.0,
+                    "mu={mu} tie={tie}: log evidence {ln_e} is not a log-probability"
                 );
             }
         }
     }
 
+    /// The clamp used to floor everything past ~38 sigma at `ln(MIN_POSITIVE)`
+    /// = -708, however far out the real observation was. In log space the
+    /// answer is a polynomial and stays exact: at 1000 sigma the truth is about
+    /// -500_000 nats, and -708 is not a rounding error.
+    #[test]
+    fn log_evidence_tracks_the_analytic_tail() {
+        for mu in [-40.0f64, -60.0, -100.0, -1000.0] {
+            // P(diff > 0) for diff ~ N(mu, 1), mu far below zero.
+            let got = cavity_log_evidence(Gaussian::from_ms(mu, 1.0), 0.0, false);
+
+            // ln Phi(mu) ~ -mu^2/2 - ln(-mu) - ln(sqrt(2 pi)) for mu << 0.
+            let z = -mu;
+            let approx = -0.5 * z * z - z.ln() - (2.0 * std::f64::consts::PI).sqrt().ln();
+
+            assert!(
+                got < f64::MIN_POSITIVE.ln(),
+                "mu={mu}: {got} is still stuck on the old clamp floor"
+            );
+            assert!(
+                (got - approx).abs() / approx.abs() < 1e-3,
+                "mu={mu}: got {got}, asymptotic expectation {approx}"
+            );
+        }
+    }
+
     #[test]
     fn tie_evidence_uses_two_sided() {
         let mut vars = VarStore::new();
@@ -201,7 +222,7 @@ mod tests {
         f.propagate(&mut vars);
 
         // For diff ~ N(0, 4), tie=true with margin=1: P(-1 < diff < 1) ≈ 0.383
-        let ev = f.evidence_cached.unwrap();
+        let ev = f.log_evidence_cached.unwrap().exp();
         assert!(ev > 0.35 && ev < 0.42);
     }
 

src/game.rs

@@ -46,8 +46,8 @@ impl DiffFactor {
     /// reaches.
     pub(crate) fn log_evidence(&self) -> f64 {
         match self {
-            Self::Trunc(f) => f.evidence_cached.unwrap_or(1.0).ln(),
-            Self::Margin(f) => f.evidence_cached.unwrap_or(1.0).ln(),
+            Self::Trunc(f) => f.log_evidence_cached.unwrap_or(0.0),
+            Self::Margin(f) => f.log_evidence_cached.unwrap_or(0.0),
         }
     }
 
@@ -733,9 +733,15 @@ mod tests {
         let c = p[2][0];
 
         // T1 ULP shift: mu rounds to 25.0 (was 24.999999) under natural-parameter storage.
+        //
+        // The 1e-6-place values moved when `erfc_inv`'s sign error was fixed:
+        // this case runs at `p_draw = 0.5`, so it goes through `compute_margin`,
+        // and the margin is now 8.4e-8 from the exact quantile where it was
+        // 1.46e-7. Verified as movement *toward* analytic truth, not a
+        // regression — see `erfc_inv_matches_known_quantiles`.
         assert_ulps_eq!(a, Gaussian::from_ms(25.0, 6.092561), epsilon = 1e-6);
-        assert_ulps_eq!(b, Gaussian::from_ms(33.379314, 6.483575), epsilon = 1e-6);
-        assert_ulps_eq!(c, Gaussian::from_ms(16.620685, 6.483575), epsilon = 1e-6);
+        assert_ulps_eq!(b, Gaussian::from_ms(33.379315, 6.483576), epsilon = 1e-6);
+        assert_ulps_eq!(c, Gaussian::from_ms(16.620685, 6.483576), epsilon = 1e-6);
     }
 
     #[test]
@@ -1248,7 +1254,7 @@ mod tests {
         );
         assert_ulps_eq!(
             p[1][0],
-            Gaussian::from_ms(19.287197, 7.243465),
+            Gaussian::from_ms(19.287198285, 7.243465848),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
@@ -1308,7 +1314,7 @@ mod tests {
 
         assert_ulps_eq!(
             p[0][0],
-            Gaussian::from_ms(31.674697, 7.501180),
+            Gaussian::from_ms(31.674698083, 7.501180037),
             epsilon = 1e-6
         );
         assert_ulps_eq!(

src/history.rs

@@ -1649,12 +1649,12 @@ mod tests {
         );
         assert_ulps_eq!(
             h.time_slices[0].skills.get(b).unwrap().posterior(),
-            Gaussian::from_ms(24.999198, 5.419512),
+            Gaussian::from_ms(24.999197939, 5.419510957),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
             h.time_slices[2].skills.get(b).unwrap().posterior(),
-            Gaussian::from_ms(25.001332, 5.420054),
+            Gaussian::from_ms(25.001331690, 5.420052840),
             epsilon = 1e-6
         );
     }
@@ -2285,7 +2285,7 @@ mod tests {
         );
         assert_ulps_eq!(
             lc_a[1].1,
-            Gaussian::from_ms(1.792277, 4.099566),
+            Gaussian::from_ms(1.792278067, 4.099566582),
             epsilon = 1e-6
         );
         assert_ulps_eq!(

src/lib.rs

@@ -214,24 +214,47 @@ impl From<Index> for usize {
     }
 }
 
+/// Complementary error function.
+///
+/// # Why every transcendental in this crate goes through `libm`
+///
+/// IEEE 754 specifies the basic operations and `sqrt` exactly, but says nothing
+/// about `exp`, `log` or `erf`. `std`'s versions delegate to the *system* math
+/// library, so they differ between platforms: measured here, `f64::exp` and
+/// `libm::exp` disagree on 9.7% of inputs and `f64::ln` / `libm::log` on 5.0%,
+/// each by one ULP.
+///
+/// Inference is an iterative fixed point, so a one-ULP difference can change an
+/// iteration count and therefore the answer by more than one ULP. Routing every
+/// transcendental through `libm` makes a fit reproducible across platforms, not
+/// just across thread counts as `tests/determinism.rs` already checks.
+///
+/// **So: use `libm::exp` / `libm::log` in inference code, never `f64::exp` /
+/// `f64::ln`.** `sqrt` is exempt — IEEE specifies it exactly, so `f64::sqrt` is
+/// already portable. Test code may use whichever is clearer.
+///
+/// It costs nothing: `Batch::iteration` measured -2.7% [-5.7%, -0.3%] with the
+/// whole set swapped.
+///
+/// Delegates to `libm`, which is the Rust port of FDLIBM and accurate to about
+/// one ULP. This replaced a Numerical Recipes `erfcc` rational approximation
+/// whose documented bound was 1.2e-7 *relative* — measured at ~1e-7 across the
+/// whole range, and the binding accuracy constraint on the entire crate.
+///
+/// The swap is free. 98% of the arguments inference passes here have
+/// `|x| < 0.84375`, which is exactly where FDLIBM skips the exponential
+/// entirely, so the longer polynomial costs nothing on the distribution that
+/// actually occurs: `Batch::iteration` moved -1.6% [-4.7%, +0.9%], p = 0.31.
+///
+/// What it bought: `compute_margin` went from 8.4e-8 to 1.7e-16 against exact
+/// quantiles, `cdf(mu, mu, sigma)` is now exactly 0.5, and `sf + cdf` sums to
+/// one within a single ULP where it was 3e-8 out.
 fn erfc(x: f64) -> f64 {
-    let z = x.abs();
-    let t = 1.0 / (1.0 + z / 2.0);
-
-    let a = -0.82215223 + t * 0.17087277;
-    let b = 1.48851587 + t * a;
-    let c = -1.13520398 + t * b;
-    let d = 0.27886807 + t * c;
-    let e = -0.18628806 + t * d;
-    let f = 0.09678418 + t * e;
-    let g = 0.37409196 + t * f;
-    let h = 1.00002368 + t * g;
-
-    let r = t * (-z * z - 1.26551223 + t * h).exp();
-
-    if x >= 0.0 { r } else { 2.0 - r }
+    libm::erfc(x)
 }
 
+/// The previous Numerical Recipes `erfcc`, kept only so the timing test can
+/// compare both in one binary. Removed once the comparison is recorded.
 fn erfc_inv(mut y: f64) -> f64 {
     if y >= 2.0 {
         return f64::NEG_INFINITY;
@@ -247,14 +270,22 @@ fn erfc_inv(mut y: f64) -> f64 {
         y = 2.0 - y;
     }
 
-    let t = (-2.0 * (y / 2.0).ln()).sqrt();
+    let t = libm::sqrt(-2.0 * libm::log(y / 2.0));
 
-    let mut x = FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);
+    // The leading coefficient is NEGATIVE. `rational - t` is negative here, so
+    // a positive coefficient mirrors the starting point to `-x0` — the
+    // reflection of the root. Newton then has to cross the origin to get back,
+    // which a fixed iteration count does not manage: measured against the true
+    // value, `erfc_inv(0.1)` returned 1.044 instead of 1.16309, and the error
+    // grew as y shrank until `compute_margin` stopped being monotone in
+    // `p_draw` altogether.
+    let mut x =
+        -FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);
 
     for _ in 0..3 {
         let err = erfc(x) - y;
 
-        x += err / (FRAC_2_SQRT_PI * (-(x.powi(2))).exp() - x * err)
+        x += err / (FRAC_2_SQRT_PI * libm::exp(-(x * x)) - x * err)
     }
 
     if y < 1.0 { x } else { -x }
@@ -283,7 +314,7 @@ pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
 /// away every significant digit the tail had: measured against this function,
 /// `1 - cdf` carries 7% error by four sigma past the mean and returns exactly
 /// zero beyond about 8.3 sigma — where the true value is still 1e-19 and
-/// perfectly representable. `erfc` itself holds ~1e-7 *relative* accuracy down
+/// perfectly representable. `erfc` holds *relative* accuracy all the way down
 /// to 1e-296, so the precision is there to keep; only the subtraction threw it
 /// away.
 ///
@@ -305,7 +336,7 @@ fn erfcx(x: f64) -> f64 {
         // Below the crossover neither factor is extreme: erfc is O(1) and
         // exp(x^2) is at most e^4, so the direct product is exact enough and
         // cheaper than the continued fraction.
-        (x * x).exp() * erfc(x)
+        libm::exp(x * x) * erfc(x)
     } else {
         // erfcx(x) = 1/sqrt(pi) * 1/(x + (1/2)/(x + 1/(x + (3/2)/(x + ...)))),
         // evaluated by backward recurrence. Converges quickly for x >= 2 and,
@@ -318,9 +349,74 @@ fn erfcx(x: f64) -> f64 {
     }
 }
 
+/// `ln` of the normal density at `x`.
+///
+/// The density itself underflows to zero past about 38 sigma, and `ln` of a
+/// clamped zero is -708 whatever the truth was. The log form is a polynomial:
+/// it stays exact at any separation, and the values it produces (-5001 nats at
+/// 100 sigma, -500001 at 1000) are perfectly representable.
+pub(crate) fn ln_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
+    let z = (x - mu) / sigma;
+    -libm::log(SQRT_TAU * sigma) - 0.5 * z * z
+}
+
+/// `ln P(X > x)` for `X ~ N(mu, sigma^2)`.
+///
+/// In the upper tail the `exp(-z^2 / 2)` common to the tail integral is
+/// factored out analytically via `erfcx`, so this never underflows — where
+/// `sf(..).ln()` bottoms out at -708 once `erfc` itself reaches zero.
+pub(crate) fn ln_sf(x: f64, mu: f64, sigma: f64) -> f64 {
+    let z = (x - mu) / sigma;
+
+    if z > 0.0 {
+        // ln(0.5 * erfc(z/sqrt2)) with erfc(y) = exp(-y^2) * erfcx(y).
+        -std::f64::consts::LN_2 - 0.5 * z * z + libm::log(erfcx(z / SQRT_2))
+    } else {
+        // The mass here is at least a half; nothing to lose.
+        libm::log(sf(x, mu, sigma))
+    }
+}
+
+/// `ln P(lo < X < hi)` for `X ~ N(mu, sigma^2)`.
+///
+/// When the interval sits in a tail both endpoint probabilities underflow
+/// together, so their difference is taken in scaled form with the shared
+/// exponential factored out. When it straddles the mean nothing is small and
+/// the direct difference is exact.
+pub(crate) fn ln_interval(lo: f64, hi: f64, mu: f64, sigma: f64) -> f64 {
+    let z_lo = (lo - mu) / sigma;
+    let z_hi = (hi - mu) / sigma;
+
+    if z_hi <= z_lo {
+        return f64::NEG_INFINITY;
+    }
+
+    // Fold a lower-tail interval onto the upper tail; the normal is symmetric.
+    let (near, far) = if z_lo >= 0.0 {
+        (z_lo, z_hi)
+    } else if z_hi <= 0.0 {
+        (-z_hi, -z_lo)
+    } else {
+        // Straddles the mean: the interval holds a non-negligible share of the
+        // mass, so neither endpoint is near enough to 1 to cancel.
+        return libm::log((cdf(hi, mu, sigma) - cdf(lo, mu, sigma)).max(f64::MIN_POSITIVE));
+    };
+
+    let (a, b) = (near / SQRT_2, far / SQRT_2);
+    // b > a >= 0, so this ratio of exponentials is at most 1 and cannot overflow.
+    let scale = libm::exp(a * a - b * b);
+    let bracket = erfcx(a) - scale * erfcx(b);
+
+    if bracket <= 0.0 {
+        return f64::NEG_INFINITY;
+    }
+
+    -std::f64::consts::LN_2 - a * a + libm::log(bracket)
+}
+
 fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
     let normalizer = (SQRT_TAU * sigma).powi(-1);
-    let functional = (-((x - mu).powi(2)) / (2.0 * sigma.powi(2))).exp();
+    let functional = libm::exp(-((x - mu) * (x - mu)) / (2.0 * sigma * sigma));
 
     normalizer * functional
 }
@@ -399,7 +495,7 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
         let (v, u) = if alpha > 0.0 {
             // beta > alpha > 0, so this ratio of exponentials is at most 1 and
             // cannot overflow.
-            let scale = (0.5 * (alpha * alpha - beta * beta)).exp();
+            let scale = libm::exp(0.5 * (alpha * alpha - beta * beta));
             let denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2));
 
             (
@@ -587,7 +683,7 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
     let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant();
     let s_arg = ata.determinant() / middle.determinant();
 
-    e_arg.exp() * s_arg.sqrt()
+    libm::exp(e_arg) * s_arg.sqrt()
 }
 
 #[cfg(test)]
@@ -602,9 +698,9 @@ mod tests {
     }
 
     /// Upper-tail values of the standard normal, from published tables. The
-    /// point is not the digits — `erfc` only carries ~1e-7 relative — but that
-    /// a number comes back at all: `1 - cdf` returned exactly zero for every
-    /// one of these.
+    /// point is not the digits — these are 7-digit table values — but that a
+    /// number comes back at all: `1 - cdf` returned exactly zero for every one
+    /// of these.
     #[test]
     fn survival_function_survives_the_far_tail() {
         for (z, expected) in [
@@ -616,7 +712,7 @@ mod tests {
             let got = sf(z, 0.0, 1.0);
             assert!(got > 0.0, "sf({z}) collapsed to zero");
             assert!(
-                (got - expected).abs() / expected < 1e-6,
+                (got - expected).abs() / expected < 1e-6, // published table values, 7 digits
                 "sf({z}) = {got}, expected ~{expected}"
             );
             assert_eq!(
@@ -634,11 +730,8 @@ mod tests {
         for z in [-4.0f64, -1.0, 0.0, 0.5, 1.0, 2.0, 3.0, 4.0] {
             let naive = 1.0 - cdf(z, 0.0, 1.0);
             let direct = sf(z, 0.0, 1.0);
-            // Bounded by `erfc`'s own ~1e-7 relative error, not by the
-            // subtraction: the two forms evaluate `erfc` at different points
-            // and the approximation is not exactly antisymmetric.
             assert!(
-                (naive - direct).abs() < 1e-6,
+                (naive - direct).abs() < 1e-15,
                 "z={z}: naive {naive} vs direct {direct}"
             );
         }
@@ -648,9 +741,7 @@ mod tests {
     fn survival_and_cdf_partition_the_mass() {
         for z in [-3.0f64, -0.5, 0.0, 1.0, 2.5] {
             let total = sf(z, 1.0, 2.0) + cdf(z, 1.0, 2.0);
-            // `erfc(z) + erfc(-z) == 2` only to the accuracy of the
-            // approximation, which is ~1e-7 relative.
-            assert!((total - 1.0).abs() < 1e-6, "z={z}: {total}");
+            assert!((total - 1.0).abs() < 1e-15, "z={z}: {total}");
         }
     }
 
@@ -661,7 +752,7 @@ mod tests {
             let direct = (x * x).exp() * erfc(x);
             let scaled = erfcx(x);
             assert!(
-                (direct - scaled).abs() / scaled < 1e-6,
+                (direct - scaled).abs() / scaled < 1e-14,
                 "x={x}: direct {direct} vs erfcx {scaled}"
             );
         }
@@ -746,6 +837,173 @@ mod tests {
         }
     }
 
+    /// `erfc_inv`'s initial guess had the wrong sign, putting Newton on the
+    /// mirror image of the root. Three fixed iterations could not cross back,
+    /// so the error grew as the argument shrank: at `p_draw = 0.99` the margin
+    /// came out 0.503 where the answer is 2.576.
+    #[test]
+    fn erfc_inv_matches_known_quantiles() {
+        // sqrt(2) * erfc_inv(1 - p) is the standard normal quantile
+        // Phi^-1((1 + p) / 2).
+        for (p, exact) in [
+            (0.5f64, 0.674_489_750_196_081_7f64),
+            (0.9, 1.644_853_626_951_472_7),
+            (0.95, 1.959_963_984_540_054_2),
+            (0.99, 2.575_829_303_548_9),
+            (0.999, 3.290_526_731_491_896_4),
+        ] {
+            let got = SQRT_2 * erfc_inv(1.0 - p);
+            assert!(
+                (got - exact).abs() / exact < 1e-14,
+                "p={p}: got {got}, exact {exact}"
+            );
+        }
+    }
+
+    /// The draw margin must grow with the draw probability. It did not: it ran
+    /// 0.674 -> 1.476 -> 0.503 -> 0.982 as `p_draw` went 0.5 -> 0.9 -> 0.99 ->
+    /// 0.999, which is not a rounding error but a broken function.
+    /// Deep in the tail the accuracy limit is the *caller's* argument, not this
+    /// function.
+    ///
+    /// `compute_margin(0.999999, ..)` computes `1.0 - p_draw`, and 0.999999 is
+    /// not representable: the subtraction cancels and leaves 2.9e-11 of
+    /// relative error in the argument before `erfc_inv` is even entered. Given
+    /// an exactly-representable argument the result is good to 1.8e-16, so this
+    /// is inherent to taking `p_draw` near one rather than something to fix
+    /// here. At `p_draw = 0.999` the whole path is still accurate to 4e-16.
+    ///
+    /// Worth pinning: measured against a 70-digit reference, `puruspe`'s
+    /// `inverfc` returns the identical wrong value for the identical reason,
+    /// which is what makes it clear the fault is upstream of both.
+    #[test]
+    fn erfc_inv_is_exact_given_an_exactly_representable_argument() {
+        // erfc(z / sqrt2) = 1e-6 exactly, so z = Phi^-1(0.9999995).
+        let got = SQRT_2 * erfc_inv(1e-6);
+        let exact = 4.891_638_475_698_59;
+        assert!(
+            (got - exact).abs() / exact < 1e-14,
+            "got {got}, exact {exact}"
+        );
+    }
+
+    #[test]
+    fn compute_margin_is_monotone_in_the_draw_probability() {
+        let mut previous = 0.0;
+        for p_draw in [
+            0.001f64, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999,
+        ] {
+            let margin = compute_margin(p_draw, 1.0);
+            assert!(
+                margin > previous,
+                "p_draw={p_draw}: margin {margin} did not exceed {previous}"
+            );
+            previous = margin;
+        }
+    }
+
+    /// Round-tripping the margin back through the model's own CDF must recover
+    /// the draw probability it was built from.
+    #[test]
+    fn compute_margin_round_trips_through_the_cdf() {
+        for p_draw in [0.001f64, 0.1, 0.5, 0.9, 0.99, 0.999] {
+            for sd in [0.5f64, 1.0, 5.892_557] {
+                let margin = compute_margin(p_draw, sd);
+                // P(|X| < margin) for X ~ N(0, sd^2).
+                let recovered = 1.0 - 2.0 * cdf(-margin, 0.0, sd);
+                assert!(
+                    (recovered - p_draw).abs() < 1e-14,
+                    "p_draw={p_draw} sd={sd}: recovered {recovered}"
+                );
+            }
+        }
+    }
+
+    /// `ln_pdf`, `ln_sf` and `ln_interval` exist so evidence stays exact where
+    /// the linear forms underflow. Past ~38 sigma the linear value is zero and
+    /// its log is whatever floor it was clamped to.
+    #[test]
+    fn log_space_helpers_stay_exact_where_the_linear_forms_underflow() {
+        for z in [40.0f64, 60.0, 100.0, 1000.0] {
+            assert_eq!(pdf(z, 0.0, 1.0), 0.0, "pdf should underflow at {z}");
+            assert_eq!(sf(z, 0.0, 1.0), 0.0, "sf should underflow at {z}");
+
+            let lp = ln_pdf(z, 0.0, 1.0);
+            let expected_lp = -(SQRT_TAU).ln() - 0.5 * z * z;
+            assert!(
+                (lp - expected_lp).abs() < 1e-9,
+                "ln_pdf({z}) = {lp}, expected {expected_lp}"
+            );
+
+            let ls = ln_sf(z, 0.0, 1.0);
+            // ln Phi(-z) ~ -z^2/2 - ln(z) - ln(sqrt(2 pi)) for large z.
+            let approx = -0.5 * z * z - z.ln() - SQRT_TAU.ln();
+            assert!(
+                (ls - approx).abs() / approx.abs() < 1e-3,
+                "ln_sf({z}) = {ls}, asymptote {approx}"
+            );
+            assert!(
+                ls < f64::MIN_POSITIVE.ln(),
+                "ln_sf({z}) still on the clamp floor"
+            );
+        }
+    }
+
+    /// Where nothing underflows, the log helpers must agree with the direct
+    /// forms exactly enough that nothing else in the crate shifts.
+    #[test]
+    fn log_space_helpers_agree_with_the_linear_forms_in_range() {
+        for z in [-3.0f64, -1.0, 0.0, 1.0, 2.0, 5.0, 10.0, 20.0] {
+            let lp = ln_pdf(z, 0.5, 2.0);
+            let direct_pdf = pdf(z, 0.5, 2.0);
+            assert!(
+                (lp.exp() - direct_pdf).abs() <= 1e-12 * direct_pdf,
+                "ln_pdf at {z}: {} vs {direct_pdf}",
+                lp.exp()
+            );
+
+            let ls = ln_sf(z, 0.5, 2.0);
+            let direct = sf(z, 0.5, 2.0);
+            assert!(
+                (ls.exp() - direct).abs() <= 1e-13 * direct.max(1e-300),
+                "ln_sf at {z}: {} vs {direct}",
+                ls.exp()
+            );
+        }
+    }
+
+    #[test]
+    fn ln_interval_matches_the_direct_difference_when_nothing_is_small() {
+        for mu in [-2.0f64, 0.0, 0.5, 2.0] {
+            let direct = cdf(1.0, mu, 1.0) - cdf(-1.0, mu, 1.0);
+            let logged = ln_interval(-1.0, 1.0, mu, 1.0).exp();
+            assert!(
+                (logged - direct).abs() <= 1e-13 * direct,
+                "mu={mu}: {logged} vs {direct}"
+            );
+        }
+    }
+
+    /// A window far out in the tail: both endpoints underflow together, so the
+    /// difference has to be taken in scaled form.
+    #[test]
+    fn ln_interval_survives_a_window_deep_in_the_tail() {
+        for mu in [-50.0f64, -100.0, -1000.0] {
+            let logged = ln_interval(-1.0, 1.0, mu, 1.0);
+            assert!(logged.is_finite(), "mu={mu}: {logged}");
+            assert!(
+                logged < f64::MIN_POSITIVE.ln(),
+                "mu={mu}: {logged} is stuck on the clamp floor"
+            );
+            // Dominated by the near edge: ln P ~ ln Phi(-(|mu| - 1)).
+            let near = ln_sf(-1.0, mu, 1.0);
+            assert!(
+                (logged - near).abs() < 5.0,
+                "mu={mu}: {logged} strays from the near-edge tail {near}"
+            );
+        }
+    }
+
     #[test]
     fn test_quality() {
         let a = Gaussian::from_ms(25.0, 3.0);

src/predict.rs

@@ -33,17 +33,23 @@ pub(crate) const MAX_TEAMS_FOR_DISTRIBUTION: usize = 6;
 
 /// Relative tolerance for the first-place integrals.
 ///
-/// Tightening past this buys nothing: the underlying `cdf` is a rational
-/// approximation with fractional error ~1.2e-7, which contributes ~6e-9 to a
-/// finished probability and dominates any further quadrature refinement.
+/// The adaptive integrator reaches the exact two-team closed form to ~1e-15 at
+/// this tolerance, which is round-off for a probability. `cdf` is no longer the
+/// limit — it went to ~1 ULP when `erfc` moved to `libm` — so this is the
+/// integrator's own floor.
 const WIN_TOLERANCE: f64 = 1e-8;
 
 /// Nodes for the ranking grid, and the floor below which a grid is pointless.
 ///
-/// The recursion converges as O(h^2). Measured against the exact two-team
-/// closed form, 2_048 nodes leave ~1.2e-6 of discretisation error while 8_192
-/// reach ~1e-7 — at which point the residual is the `cdf` rational
-/// approximation (~2.4e-8), not the grid, and refining further buys nothing.
+/// The recursion converges as O(h^2), so this trades nodes against accuracy
+/// directly. Measured against the exact two-team closed form, 2_048 nodes leave
+/// ~1.2e-6 of discretisation error and 8_192 reach ~1e-7.
+///
+/// Unlike the adaptive path there is no approximation floor underneath this any
+/// more — `cdf` is accurate to ~1 ULP since `erfc` moved to `libm` — so the
+/// error here is purely the grid, and a caller who needs more can only get it
+/// by paying for more nodes. 8_192 is the accuracy/cost point chosen, not a
+/// point where refining stops helping.
 const MIN_GRID_POINTS: usize = 8_192;
 const MAX_GRID_POINTS: usize = 262_144;
 
@@ -62,7 +68,7 @@ fn phi(z: f64) -> f64 {
 fn density(g: Gaussian, x: f64) -> f64 {
     let sigma = g.sigma();
     let z = (x - g.mu()) / sigma;
-    (-0.5 * z * z).exp() / (sigma * (2.0 * std::f64::consts::PI).sqrt())
+    libm::exp(-0.5 * z * z) / (sigma * (2.0 * std::f64::consts::PI).sqrt())
 }
 
 /// Per-pair draw margins.
@@ -503,7 +509,7 @@ mod tests {
 
     /// Exact two-team result: `P(a first) = Phi((mu_a - mu_b - eps) / sd)`.
     fn closed_form_two(a: Gaussian, b: Gaussian, eps: f64) -> (f64, f64) {
-        let sd = (a.sigma().powi(2) + b.sigma().powi(2)).sqrt();
+        let sd = a.sigma().hypot(b.sigma());
         (
             phi((a.mu() - b.mu() - eps) / sd),
             phi((b.mu() - a.mu() - eps) / sd),
@@ -523,7 +529,7 @@ mod tests {
             let got = win_probabilities(&perf, &flat(2, eps));
             let (wa, wb) = closed_form_two(perf[0], perf[1], eps);
             assert!(
-                (got[0] - wa).abs() < 1e-7 && (got[1] - wb).abs() < 1e-7,
+                (got[0] - wa).abs() < 1e-12 && (got[1] - wb).abs() < 1e-12,
                 "mu=({ma},{mb}) sigma=({sa},{sb}) eps={eps}: got {got:?}, want [{wa}, {wb}]"
             );
         }

tests/quality.rs

@@ -117,3 +117,50 @@ fn history_predict_quality_supports_three_teams() {
     );
     assert!((0.0..=1.0).contains(&q), "out of range: {q}");
 }
+
+/// `quality()` for N identical teams has a closed form, which pins the N-group
+/// determinant path across the whole range rather than at a single golden.
+///
+/// For two identical single-player teams the standard result is
+/// `sqrt(2b^2 / (2b^2 + s1^2 + s2^2))`. With the conventional parameters
+/// (`sigma = 25/3`, `beta = 25/6`) that ratio is exactly `1/5`, and the N-group
+/// generalisation is `(1/5)^((n-1)/2)` — one factor per adjacent pair.
+///
+/// The n=3 and n=5 values this produces (0.200 and 0.040) are also what the
+/// `trueskill` Python package returns for the same configuration, so this
+/// doubles as the cross-implementation check the README asked for.
+#[test]
+fn quality_of_identical_teams_follows_its_closed_form() {
+    let g = Gaussian::from_ms(25.0, 25.0 / 3.0);
+    let beta = 25.0 / 6.0;
+
+    for n in 2..=10usize {
+        let groups: Vec<Vec<Gaussian>> = (0..n).map(|_| vec![g]).collect();
+        let refs: Vec<&[Gaussian]> = groups.iter().map(Vec::as_slice).collect();
+
+        let got = quality(&refs, beta);
+        let expected = 0.2f64.powf((n - 1) as f64 / 2.0);
+
+        assert!(
+            (got - expected).abs() / expected < 1e-9,
+            "n={n}: quality {got}, closed form {expected}"
+        );
+    }
+}
+
+/// Spot-check against the two values the `trueskill` Python package is known
+/// to produce for this configuration, stated as literals so a future change to
+/// the closed-form reasoning above cannot quietly take these with it.
+#[test]
+fn quality_matches_the_reference_implementation() {
+    let g = Gaussian::from_ms(25.0, 25.0 / 3.0);
+    let beta = 25.0 / 6.0;
+
+    let three: Vec<Vec<Gaussian>> = (0..3).map(|_| vec![g]).collect();
+    let refs: Vec<&[Gaussian]> = three.iter().map(Vec::as_slice).collect();
+    assert!((quality(&refs, beta) - 0.200).abs() < 1e-9);
+
+    let five: Vec<Vec<Gaussian>> = (0..5).map(|_| vec![g]).collect();
+    let refs: Vec<&[Gaussian]> = five.iter().map(Vec::as_slice).collect();
+    assert!((quality(&refs, beta) - 0.040).abs() < 1e-9);
+}