fix: route every transcendental through libm, and combine sigmas with hypot

Follow-on from #41, which added `libm` for `erfc`. Surveying what else
the dependency offers: its unique surface over `std` is `erf`/`erfc`,
`lgamma`/`tgamma` and Bessel functions, and only the first was ever
needed. But the survey found something better than another special
function.

`std`'s `exp` and `ln` delegate to the *system* math library. IEEE 754
specifies the basic operations and `sqrt` exactly and says nothing about
transcendentals, so those differ per platform. Measured here over 200k
inputs:

    exp: 19425/200000 differ from libm (worst 1 ulp)
    log:  9932/200000 differ

Inference is an iterative fixed point, so a one-ULP difference can change
an iteration count and move the answer by more than one ULP. Routing
every transcendental through `libm` makes a fit reproducible across
platforms — a stronger guarantee than `tests/determinism.rs`, which only
covers thread counts.

It costs nothing. `Batch::iteration` measured -2.7% [-5.7%, -0.3%] with
the whole set swapped, and not one golden moved.

Also switches the two places that combined sigmas as
`sqrt(a^2 + b^2)` to `hypot`. Squaring overflows to infinity above
~1.3e154 and flushes to zero below ~1.5e-154 — measured, the naive form
returns `inf` where `hypot` returns 1.41e160 — and `Gaussian`'s
constructors are public, so a caller can reach both ends.

Deliberately not done: rewriting the KL divergence's `ln` of a ratio via
`ln_1p`. The cancellation is real as the ratio approaches one, but
measured absolute error is at most ~1e-11 in a quantity of order 0.4
nats, so it changes nothing.

The invariant is recorded in `CLAUDE.md` and on `erfc`'s own docs, since
nothing enforces it mechanically.

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

src/lib.rs

@@ -216,6 +216,26 @@ 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
@@ -250,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
@@ -265,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 }
@@ -316,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,
@@ -337,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)`.
@@ -350,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))
     }
 }
 
@@ -379,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))
-            .max(f64::MIN_POSITIVE)
-            .ln();
+        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 = (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
 }
@@ -477,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));
 
             (
@@ -665,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)]