fix: stop destroying tail precision in evidence and truncation

`erfc` is sound — it holds ~1e-7 *relative* accuracy down to 1e-296 with
no tail degradation. Three expressions built on it threw that away by
subtracting quantities that both approach the same value.

1. `cavity_evidence` computed `1.0 - cdf(margin, ..)`, which is
   algebraically `sf(margin, ..)` and numerically a catastrophe: 7%
   error by eight sigma, and exactly zero past ~8.3, where the true
   probability is 1e-19 and perfectly representable. Clamped, that
   reached `log_evidence` as ln(f64::MIN_POSITIVE) = -708 whatever the
   truth was — off by 665 nats at nine sigma.

   `1 - cdf` is smallest precisely when the result contradicts the
   prior, so the model-comparison number was worst for upsets: the
   observation it exists to notice. Adds `sf`, the survival function,
   computed without the subtraction. The tie branch picks whichever tail
   keeps both of its terms small, for the same reason.

2. `v_w` computed the inverse Mills ratio as `pdf(-a) / cdf(-a)`. Both
   underflow together past about 39 sigma, giving `0 / 0` and putting
   NaN straight into the posterior. Adds `erfcx`, so the shared
   `exp(-alpha^2 / 2)` cancels analytically instead of being evaluated
   twice and divided.

3. With that fixed, `w = v * (v - alpha)` became the next casualty: `v`
   tends to `alpha`, so the gap lost every digit and drove `w` above 1,
   making `sqrt(1 - w)` NaN at alpha = 1e6. The gap now comes from its
   asymptotic series, which forms no difference at all. The tie branch
   had the same defect one expression over — `v * v - u` with both terms
   at 1e18 returned w = -128 — and a far-tail window is
   indistinguishable from a half-line, so it shares the asymptotic.

No public signature changes, and no existing golden moved: every one of
these only alters regions the old code got wrong. The two identity tests
are asserted at 1e-6 rather than tighter because `erfc` is not exactly
antisymmetric — `erfc(z) + erfc(-z)` differs from 2 by ~3e-8, and
`erfc(0)` returns 1.00000003. That floor is tracked separately.

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

src/factor/trunc.rs

@@ -2,6 +2,7 @@ use crate::{
     N_INF, approx, cdf,
     factor::{Factor, VarId, VarStore},
     gaussian::Gaussian,
+    sf,
 };
 
 /// EP truncation factor on a diff variable.
@@ -74,16 +75,29 @@ impl Factor for TruncFactor {
 
 /// P(diff > margin) for non-tie, P(|diff| < margin) for tie.
 ///
-/// Clamped to a positive floor: for a near-certain outcome the tail rounds to
+/// Both branches pick whichever tail keeps their terms *small*, because the
-/// exactly 0.0, and the `erfc` approximation used by `cdf` carries ~1e-7 error
+/// alternative is subtracting two numbers that both approach 1. That
-/// so it can even return slightly more than 1.0, making the difference
+/// subtraction is not a rounding detail: it loses every digit of an unlikely
-/// negative. Either would send `log_evidence` to `-inf` or NaN and poison the
+/// outcome's evidence, and an unlikely outcome is precisely the one worth
-/// sum across the whole history.
+/// 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 {
+    let (mu, sigma) = (diff.mu(), diff.sigma());
+
     let raw = if tie {
-        cdf(margin, diff.mu(), diff.sigma()) - cdf(-margin, diff.mu(), diff.sigma())
+        if mu < -margin {
+            // Both CDFs sit against 1 here; both survival terms are small.
+            sf(-margin, mu, sigma) - sf(margin, mu, sigma)
+        } else {
+            cdf(margin, mu, sigma) - cdf(-margin, mu, sigma)
+        }
     } else {
-        1.0 - cdf(margin, diff.mu(), diff.sigma())
+        sf(margin, mu, sigma)
     };
 
     raw.clamp(f64::MIN_POSITIVE, 1.0)
@@ -132,6 +146,52 @@ mod tests {
         assert_eq!(f.evidence_cached.unwrap(), first);
     }
 
+    /// The defect this guards: `1 - cdf` collapsed to zero for a surprising
+    /// result, the clamp turned that into `f64::MIN_POSITIVE`, and
+    /// `log_evidence` reported ln of *that* — about -708 whatever the truth
+    /// was. An upset is the observation a model-comparison score exists to
+    /// notice, so it was wrong exactly where it mattered.
+    #[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);
+
+        assert!(
+            evidence > f64::MIN_POSITIVE,
+            "evidence collapsed onto the clamp floor: {evidence}"
+        );
+        // P(X > 0) for X ~ N(-9, 1) is the standard normal tail at 9 sigma.
+        assert!(
+            (evidence - 1.128_588e-19).abs() / 1.128_588e-19 < 1e-6,
+            "expected ~1.13e-19, got {evidence}"
+        );
+        assert!(
+            (evidence.ln() + 43.628).abs() < 1e-2,
+            "log evidence {} should be about -43.6, not -708",
+            evidence.ln()
+        );
+    }
+
+    /// 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.
+    #[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);
+                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"
+                );
+            }
+        }
+    }
+
     #[test]
     fn tie_evidence_uses_two_sided() {
         let mut vars = VarStore::new();

src/lib.rs

@@ -168,6 +168,21 @@ pub const ITERATIONS: usize = 30;
 pub const MAX_PREDICTED_TEAMS: usize = predict::MAX_TEAMS_FOR_DISTRIBUTION;
 
 const SQRT_TAU: f64 = 2.5066282746310002;
+/// `1 / sqrt(pi)`, the leading factor of the `erfcx` continued fraction.
+const FRAC_1_SQRT_PI: f64 = 0.564_189_583_547_756_3;
+/// `sqrt(2 / pi)`, the numerator of the inverse Mills ratio in scaled form.
+const SQRT_2_OVER_PI: f64 = 0.797_884_560_802_865_4;
+/// How many window widths into the tail before a tie window is treated as a
+/// half-line. Beyond this the truncated mass is concentrated within `1/alpha`
+/// of the near edge, so the far edge contributes nothing measurable.
+const HALF_LINE_WINDOW: f64 = 10.0;
+/// Where `v - alpha` switches from subtraction to its asymptotic series.
+///
+/// The subtraction loses roughly `eps * alpha^2` of relative precision, and the
+/// 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.
+const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0;
 
 pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0);
 pub const N00: Gaussian = Gaussian::from_ms(0.0, 0.0);
@@ -259,6 +274,50 @@ pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
     0.5 * erfc(z)
 }
 
+/// `P(X > x)` for `X ~ N(mu, sigma^2)`.
+///
+/// The survival function, computed directly rather than as `1 - cdf(..)`.
+///
+/// The two are algebraically identical and numerically are not. `cdf` returns
+/// a value approaching 1 for an upper tail, so subtracting it from 1 cancels
+/// 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
+/// to 1e-296, so the precision is there to keep; only the subtraction threw it
+/// away.
+///
+/// This matters most where evidence is smallest, which is exactly where an
+/// upset makes it interesting: `ln` of a clamped zero is -708 regardless of
+/// whether the truth was -43 or -600.
+pub(crate) fn sf(x: f64, mu: f64, sigma: f64) -> f64 {
+    0.5 * erfc((x - mu) / (sigma * SQRT_2))
+}
+
+/// `e^(x^2) * erfc(x)`, the scaled complementary error function, for `x >= 0`.
+///
+/// Exists so the exponential factor common to a Gaussian density and its tail
+/// integral can be cancelled *analytically* instead of being computed twice
+/// and divided. Both underflow to zero past about 26 sigma, and their ratio is
+/// then `0/0` — finite in the limit, `NaN` in floating point.
+fn erfcx(x: f64) -> f64 {
+    if x < 2.0 {
+        // 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)
+    } 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,
+        // unlike the product form, never touches an exponential.
+        let mut f = 0.0;
+        for n in (1..=60u32).rev() {
+            f = (f64::from(n) * 0.5) / (x + f);
+        }
+        FRAC_1_SQRT_PI / (x + f)
+    }
+}
+
 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();
@@ -266,25 +325,100 @@ fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
     normalizer * functional
 }
 
+/// Truncated-Gaussian correction terms `(v, w)`.
+///
+/// `v` shifts the mean and `w` shrinks the variance. Both are ratios whose
+/// numerator and denominator underflow together in the tails, so both are
+/// computed in scaled form there: the shared `exp(-alpha^2 / 2)` is cancelled
+/// analytically rather than evaluated and divided out. Without that, a
+/// truncation point beyond about 39 sigma produced `0 / 0` and put `NaN`
+/// straight into the posterior.
+/// Truncation terms for a boundary `alpha` standard deviations into the upper
+/// tail, from the asymptotic expansion of the inverse Mills ratio.
+///
+/// `v` tends to `alpha` out here, so the gap between them cannot be obtained by
+/// subtracting one from the other — the series computes the gap directly, and
+/// `w = v * gap` then never forms the difference of two large near-equal
+/// numbers. A far-tail *window* behaves like a half-line once it is more than a
+/// few multiples of its own width from the mean, so the tie branch shares this.
+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 v = alpha + gap;
+
+    (v, v * gap)
+}
+
 fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
     if !tie {
         let alpha = (margin - mu) / sigma;
 
-        let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
+        // v is the inverse Mills ratio, phi(alpha) / Phi(-alpha), and w needs
-        let w = v * (v + (-alpha));
+        // the gap `v - alpha` as well as v itself. Far into the tail v tends to
+        // alpha, so that gap is a subtraction of two nearly equal numbers and
+        // loses every digit it has: at alpha = 1e6 it drove w above 1 and made
+        // `sqrt(1 - w)` NaN. Past the crossover the gap comes from its
+        // asymptotic series instead, which has no subtraction in it.
+        if alpha >= ASYMPTOTIC_MILLS_ALPHA {
+            return half_line_truncation(alpha);
+        }
 
-        (v, w)
+        let (v, gap) = if alpha > 0.0 {
+            // Both terms carry exp(-alpha^2 / 2); in scaled form it cancels
+            // and the result stays exact however far into the tail alpha sits.
+            let v = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2);
+            (v, v - alpha)
+        } else {
+            // Phi(-alpha) >= 1/2 here, so the direct ratio loses nothing.
+            let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
+            (v, v - alpha)
+        };
+
+        (v, 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.
+        let flipped = mu > 0.0;
+        let mu = if flipped { -mu } else { mu };
+
         let alpha = (-margin - mu) / sigma;
         let beta = (margin - mu) / sigma;
 
-        let v = (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0))
+        // `w` comes out of `v * v - u`, and both terms grow as alpha^2 while
-            / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
+        // their difference stays O(1) — at alpha = 1e9 that subtraction had no
-        let u = (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0))
+        // digits left and returned w = -128, making `sqrt(1 - w)` nonsense.
-            / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
+        // 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 (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 denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2));
+
+            (
+                (1.0 - scale) / SQRT_TAU / denominator,
+                (alpha - beta * scale) / SQRT_TAU / denominator,
+            )
+        } else {
+            // The interval straddles the mean, so nothing here is small.
+            let denominator = cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0);
+
+            (
+                (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) / denominator,
+                (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) / denominator,
+            )
+        };
+
         let w = -(u - v.powi(2));
 
-        (v, w)
+        (if flipped { -v } else { v }, w)
     }
 }
 
@@ -467,6 +601,151 @@ mod tests {
         assert_eq!(sort_time(&[0i64, 1, 2, 0], true), vec![2, 1, 0, 3]);
     }
 
+    /// 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.
+    #[test]
+    fn survival_function_survives_the_far_tail() {
+        for (z, expected) in [
+            (9.0f64, 1.128_588e-19),
+            (12.0, 1.776_482e-33),
+            (20.0, 2.753_624e-89),
+            (37.0, 5.725_571e-300),
+        ] {
+            let got = sf(z, 0.0, 1.0);
+            assert!(got > 0.0, "sf({z}) collapsed to zero");
+            assert!(
+                (got - expected).abs() / expected < 1e-6,
+                "sf({z}) = {got}, expected ~{expected}"
+            );
+            assert_eq!(
+                1.0 - cdf(z, 0.0, 1.0),
+                0.0,
+                "the naive form should still be zero here"
+            );
+        }
+    }
+
+    /// Where no cancellation happens the two forms must agree exactly enough
+    /// that nothing else in the crate shifts.
+    #[test]
+    fn survival_function_matches_the_naive_form_where_that_form_works() {
+        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,
+                "z={z}: naive {naive} vs direct {direct}"
+            );
+        }
+    }
+
+    #[test]
+    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}");
+        }
+    }
+
+    /// `erfcx` switches formulation at x = 2; the two sides must meet.
+    #[test]
+    fn erfcx_is_continuous_across_its_crossover() {
+        for x in [1.90f64, 1.99, 1.999, 2.0, 2.001, 2.01, 2.10] {
+            let direct = (x * x).exp() * erfc(x);
+            let scaled = erfcx(x);
+            assert!(
+                (direct - scaled).abs() / scaled < 1e-6,
+                "x={x}: direct {direct} vs erfcx {scaled}"
+            );
+        }
+    }
+
+    /// The whole reason `erfcx` exists: it stays finite and O(1/x) exactly
+    /// where `exp(x^2)` overflows and `erfc(x)` underflows.
+    #[test]
+    fn erfcx_stays_finite_where_its_factors_do_not() {
+        for x in [27.0f64, 50.0, 1.0e3, 1.0e8] {
+            let scaled = erfcx(x);
+            assert!(scaled.is_finite() && scaled > 0.0, "erfcx({x}) = {scaled}");
+            // Asymptotically erfcx(x) -> 1 / (x * sqrt(pi)).
+            let asymptote = 1.0 / (x * std::f64::consts::PI.sqrt());
+            assert!(
+                (scaled - asymptote).abs() / asymptote < 1e-2,
+                "erfcx({x}) = {scaled} strays from its asymptote {asymptote}"
+            );
+            assert!(
+                (x * x).exp().is_infinite(),
+                "x={x} should overflow the direct form"
+            );
+        }
+    }
+
+    /// Truncation must never produce a non-finite posterior. Before the scaled
+    /// formulation these returned NaN from `0 / 0` past about 39 sigma.
+    #[test]
+    fn truncation_stays_finite_arbitrarily_far_into_the_tail() {
+        for alpha in [0.0f64, 8.0, 38.0, 40.0, 100.0, 1.0e3, 1.0e6, 1.0e9, 1.0e15] {
+            for tie in [false, true] {
+                let (v, w) = v_w(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
+                assert!(v.is_finite(), "alpha={alpha} tie={tie}: v = {v}");
+                assert!(w.is_finite(), "alpha={alpha} tie={tie}: w = {w}");
+                // sigma_trunc = sigma * sqrt(1 - w) must stay real.
+                assert!(
+                    (0.0..=1.0).contains(&w),
+                    "alpha={alpha} tie={tie}: w = {w} leaves sqrt(1 - w) imaginary"
+                );
+
+                let (mu_t, sigma_t) = trunc(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
+                assert!(
+                    mu_t.is_finite() && sigma_t.is_finite(),
+                    "alpha={alpha} tie={tie}: trunc = ({mu_t}, {sigma_t})"
+                );
+            }
+        }
+    }
+
+    /// The Mills gap switches from subtraction to series at alpha = 100. Both
+    /// are supposed to be right there; if they disagree, the crossover is in
+    /// the wrong place.
+    #[test]
+    fn the_mills_gap_series_meets_the_scaled_form() {
+        for alpha in [50.0f64, 99.0, 100.0, 101.0, 200.0] {
+            let scaled = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2) - alpha;
+            let inv = alpha.recip();
+            let inv_sq = inv * inv;
+            let series = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
+            assert!(
+                (scaled - series).abs() / series < 1e-9,
+                "alpha={alpha}: scaled {scaled} vs series {series}"
+            );
+        }
+    }
+
+    /// Folding the tie branch to `mu <= 0` is only valid if v is odd in mu and
+    /// w is even. Assert the symmetry the implementation relies on.
+    #[test]
+    fn tie_truncation_is_odd_in_v_and_even_in_w() {
+        for mu in [0.5f64, 3.0, 20.0, 40.0, 100.0, 1.0e3] {
+            let (v_pos, w_pos) = v_w(mu, 1.0, 1.0, true);
+            let (v_neg, w_neg) = v_w(-mu, 1.0, 1.0, true);
+            assert!(
+                (v_pos + v_neg).abs() < 1e-9,
+                "mu={mu}: v should be odd, got {v_pos} and {v_neg}"
+            );
+            assert!(
+                (w_pos - w_neg).abs() < 1e-9,
+                "mu={mu}: w should be even, got {w_pos} and {w_neg}"
+            );
+        }
+    }
+
     #[test]
     fn test_quality() {
         let a = Gaussian::from_ms(25.0, 3.0);