diff --git a/CLAUDE.md b/CLAUDE.md index 4771cde..d6c3f00 100644 --- a/CLAUDE.md +++ b/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 diff --git a/src/acquisition.rs b/src/acquisition.rs index d3bc041..29ad6a2 100644 --- a/src/acquisition.rs +++ b/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. diff --git a/src/factor/margin.rs b/src/factor/margin.rs index 2d96900..dd23403 100644 --- a/src/factor/margin.rs +++ b/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 diff --git a/src/lib.rs b/src/lib.rs index 968daaa..c26971f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -216,6 +216,26 @@ impl From 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)] diff --git a/src/predict.rs b/src/predict.rs index d78ad69..2f221b3 100644 --- a/src/predict.rs +++ b/src/predict.rs @@ -68,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. @@ -509,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),