fix: keep the truncated variance representable in the far tail

`v_w` returned `w` and let `trunc` form `1 - w`. `w` tends to 1 out in
the tail, so that subtraction lost about log10(alpha^2) digits — and the
quantity it was destroying is perfectly representable.

Two separate cancellations, fixed separately.

The non-tie half: `half_line_truncation` now returns `1 - w` computed
symbolically rather than as `1 - v*gap`. With `alpha*gap = 1 - inv^2*b`
the leading ones cancel on paper instead of in floating point. Measured
against the exact truncated variance:

  alpha    before          after
  1e6      8.9e-5 rel      0.0 rel (exact)
  1e8      returns 0.0     0.0 rel (exact)

At 1e8 the old form gave `sigma_trunc = 0`, and `from_ms(mu, 0.0)` is a
point mass whose `mu()` is inf/inf = NaN. `beta(1e-8).sigma(1e-8)` with
priors 1000 apart went from Err + NaN skills to a finite fit.

The tie half is a different subtraction — `w = v^2 - u`, where both grow
as alpha^2 while their difference stays O(1). The existing escape hatch
could not cover it: it keys on `alpha * width >= HALF_LINE_WINDOW`, how
many window-widths from the mean the window sits, and a NARROW window
fails that however deep it is. Measured at alpha 1e6 with a 1e-6 window
it kept four digits and returned `1 - w = -2.4e-4` where the truth is
+2.8e-13. 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.

Over a narrow window the density is a truncated exponential in
`s = (x - alpha)/width`, whose mean and variance are closed forms, so
`v = alpha + width*m(t)` and `1 - w = width^2 * V(t)` with no large
subtraction at all. Validated against high-precision quadrature: v exact
to 4e-10, `1 - w` to 4e-10 across the region it is used in.

The crossover is on `alpha / width` rather than on either alone, because
that ratio is what says how many digits the subtraction has left — and
the approximation is most accurate exactly where the subtraction is
worst, since both improve as the window narrows.

Defaults are bit-identical (pi 0.02398318151216503 before and after).

Tests: the three reproductions from the issue, the narrow-window form
against pinned quadrature values, and a continuity sweep across all three
tie branches — a misplaced crossover is the real risk here, and a jump at
a boundary is visible even without pinning absolute values.

Closes #60

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

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 width = beta - alpha;
-            let (v, w) = half_line_truncation(alpha);
+
-            return (if flipped { -v } else { v }, w);
+        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)
 }
@@ -762,6 +849,80 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
 
 #[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

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,
@@ -163,3 +165,53 @@ fn a_nan_competitor_is_not_masked_by_a_healthy_one() {
         "{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:?}"
+        );
+    }
+}