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,12 +2,27 @@
 
 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

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.1"
+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"
 

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

@@ -74,7 +74,10 @@ impl Factor for MarginFactor {
 /// 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 {
-    let combined_sigma = (cavity.sigma().powi(2) + sigma.powi(2)).sqrt();
+    // `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);
 
     // A degenerate cavity (infinite sigma) is the only way to reach a

src/game.rs

@@ -1254,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!(
@@ -1314,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();
+    libm::erfc(x)
-    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 }
 }
 
+/// 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,7 +270,7 @@ 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));
 
     // The leading coefficient is NEGATIVE. `rational - t` is negative here, so
     // a positive coefficient mirrors the starting point to `-x0` — the
@@ -262,7 +285,7 @@ fn erfc_inv(mut y: f64) -> f64 {
     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 }
@@ -291,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.
 ///
@@ -313,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,
@@ -334,7 +357,7 @@ fn erfcx(x: f64) -> f64 {
 /// 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;
-    -(SQRT_TAU * sigma).ln() - 0.5 * z * z
+    -libm::log(SQRT_TAU * sigma) - 0.5 * z * z
 }
 
 /// `ln P(X > x)` for `X ~ N(mu, sigma^2)`.
@@ -347,10 +370,10 @@ pub(crate) fn ln_sf(x: f64, mu: f64, sigma: f64) -> f64 {
 
     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 + erfcx(z / SQRT_2).ln()
+        -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.
-        sf(x, mu, sigma).ln()
+        libm::log(sf(x, mu, sigma))
     }
 }
 
@@ -376,26 +399,24 @@ pub(crate) fn ln_interval(lo: f64, hi: f64, mu: f64, sigma: f64) -> f64 {
     } 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 (cdf(hi, mu, sigma) - cdf(lo, mu, sigma))
+        return libm::log((cdf(hi, mu, sigma) - cdf(lo, mu, sigma)).max(f64::MIN_POSITIVE));
-            .max(f64::MIN_POSITIVE)
-            .ln();
     };
 
     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 = (a * a - b * b).exp();
+    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 + bracket.ln()
+    -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
 }
@@ -474,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));
 
             (
@@ -662,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)]
@@ -677,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
+    /// point is not the digits — these are 7-digit table values — but that a
-    /// a number comes back at all: `1 - cdf` returned exactly zero for every
+    /// number comes back at all: `1 - cdf` returned exactly zero for every one
-    /// one of these.
+    /// of these.
     #[test]
     fn survival_function_survives_the_far_tail() {
         for (z, expected) in [
@@ -691,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!(
@@ -709,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}"
             );
         }
@@ -723,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
+            assert!((total - 1.0).abs() < 1e-15, "z={z}: {total}");
-            // approximation, which is ~1e-7 relative.
-            assert!((total - 1.0).abs() < 1e-6, "z={z}: {total}");
         }
     }
 
@@ -736,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}"
             );
         }
@@ -838,7 +854,7 @@ mod tests {
         ] {
             let got = SQRT_2 * erfc_inv(1.0 - p);
             assert!(
-                (got - exact).abs() / exact < 1e-6,
+                (got - exact).abs() / exact < 1e-14,
                 "p={p}: got {got}, exact {exact}"
             );
         }
@@ -847,6 +863,30 @@ mod tests {
     /// 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;
@@ -872,7 +912,7 @@ mod tests {
                 // 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-6,
+                    (recovered - p_draw).abs() < 1e-14,
                     "p_draw={p_draw} sd={sd}: recovered {recovered}"
                 );
             }
@@ -922,14 +962,10 @@ mod tests {
                 lp.exp()
             );
 
-            // Bounded by `erfc`'s ~1e-7, not tighter: above x = 2 `erfcx` uses
-            // a continued fraction accurate to ~1e-15, so the log path is the
-            // *more* accurate of the two and they part company at `erfc`'s
-            // error rather than at round-off.
             let ls = ln_sf(z, 0.5, 2.0);
             let direct = sf(z, 0.5, 2.0);
             assert!(
-                (ls.exp() - direct).abs() <= 1e-6 * direct.max(1e-300),
+                (ls.exp() - direct).abs() <= 1e-13 * direct.max(1e-300),
                 "ln_sf at {z}: {} vs {direct}",
                 ls.exp()
             );
@@ -941,11 +977,8 @@ mod tests {
         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();
-            // See `log_space_helpers_agree_with_the_linear_forms_in_range`:
-            // the gap here is `erfc`'s own error, and the log path is the more
-            // accurate side of it.
             assert!(
-                (logged - direct).abs() <= 1e-6 * direct,
+                (logged - direct).abs() <= 1e-13 * direct,
                 "mu={mu}: {logged} vs {direct}"
             );
         }

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
+/// The adaptive integrator reaches the exact two-team closed form to ~1e-15 at
-/// approximation with fractional error ~1.2e-7, which contributes ~6e-9 to a
+/// this tolerance, which is round-off for a probability. `cdf` is no longer the
-/// finished probability and dominates any further quadrature refinement.
+/// 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
+/// The recursion converges as O(h^2), so this trades nodes against accuracy
-/// closed form, 2_048 nodes leave ~1.2e-6 of discretisation error while 8_192
+/// directly. Measured against the exact two-team closed form, 2_048 nodes leave
-/// reach ~1e-7 — at which point the residual is the `cdf` rational
+/// ~1.2e-6 of discretisation error and 8_192 reach ~1e-7.
-/// approximation (~2.4e-8), not the grid, and refining further buys nothing.
+///
+/// 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}]"
             );
         }