From bdf33ac9da0c18ff0c7485595f48cb3d0263b5cc Mon Sep 17 00:00:00 2001 From: Dundy Pasupuleti Date: Sat, 31 Oct 2026 04:25:48 -0400 Subject: [PATCH 1/3] fix: keep Cauchy::inverse_cdf accurate in the lower tail Signed-off-by: Dundy Pasupuleti --- src/distribution/cauchy.rs | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/src/distribution/cauchy.rs b/src/distribution/cauchy.rs index 0d978b6f..e4352fa5 100644 --- a/src/distribution/cauchy.rs +++ b/src/distribution/cauchy.rs @@ -166,7 +166,16 @@ impl ContinuousCDF for Cauchy { if !(0.0..=1.0).contains(&x) { panic!("x must be in [0, 1]"); } else { - self.location + self.scale * (f64_consts::PI * (x - 0.5)).tan() + // tan(pi (x - 1/2)) = -1 / tan(pi x); using x (or 1 - x, which + // is exact for x >= 1/2) avoids the rounding of x - 1/2 that made + // every lower quantile below 2^-55 collapse to one value (gh #471) + if x == 0.5 { + self.location + } else if x < 0.5 { + self.location - self.scale / (f64_consts::PI * x).tan() + } else { + self.location + self.scale / (f64_consts::PI * (1.0 - x)).tan() + } } } } @@ -526,4 +535,19 @@ mod tests { density_util::check_continuous_distribution(&create_ok(-1.2, 3.4), -1500.0, 1500.0); density_util::check_continuous_distribution(&create_ok(-4.5, 6.7), -5000.0, 5000.0); } + + #[test] + fn test_inverse_cdf_lower_tail() { + // gh #471: the lower tail used to saturate at -1.63e16 below p = 2^-55 + let d = create_ok(0.0, 1.0); + for p in [1e-15f64, 2f64.powi(-55), 1e-100, 1e-300] { + approx::assert_relative_eq!(d.inverse_cdf(p), -1.0 / (core::f64::consts::PI * p), max_relative = 1e-14); + } + assert_ne!(d.inverse_cdf(2f64.powi(-56)), d.inverse_cdf(2f64.powi(-60))); + assert_eq!(d.inverse_cdf(0.0), f64::NEG_INFINITY); + assert_eq!(d.inverse_cdf(1.0), f64::INFINITY); + assert_eq!(d.inverse_cdf(0.5), 0.0); + approx::assert_relative_eq!(d.inverse_cdf(1e-3), -318.30883898555044, max_relative = 1e-14); + approx::assert_relative_eq!(d.inverse_cdf(1.0 - 1e-3), 318.30883898555044, max_relative = 1e-12); + } } From 8b27af07eec730cb5f165a25b91ebb7bf7d33932 Mon Sep 17 00:00:00 2001 From: Dundy Pasupuleti Date: Sat, 31 Oct 2026 04:45:11 -0400 Subject: [PATCH 2/3] cauchy: keep inverse_cdf exact for subnormal p below p = 1e-9, tan(pi p) equals pi p to double precision, so divide by p directly; forming pi * p rounded badly once p was subnormal. Signed-off-by: Dundy Pasupuleti --- src/distribution/cauchy.rs | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/distribution/cauchy.rs b/src/distribution/cauchy.rs index e4352fa5..ef433e2c 100644 --- a/src/distribution/cauchy.rs +++ b/src/distribution/cauchy.rs @@ -171,6 +171,10 @@ impl ContinuousCDF for Cauchy { // every lower quantile below 2^-55 collapse to one value (gh #471) if x == 0.5 { self.location + } else if x < 1e-9 { + // tan(pi x) = pi x to double precision here; dividing by x + // directly keeps a subnormal x from rounding pi x + self.location - self.scale * f64_consts::FRAC_1_PI / x } else if x < 0.5 { self.location - self.scale / (f64_consts::PI * x).tan() } else { @@ -550,4 +554,14 @@ mod tests { approx::assert_relative_eq!(d.inverse_cdf(1e-3), -318.30883898555044, max_relative = 1e-14); approx::assert_relative_eq!(d.inverse_cdf(1.0 - 1e-3), 318.30883898555044, max_relative = 1e-12); } + + #[test] + fn test_inverse_cdf_subnormal_p() { + let d = create_ok(0.0, 1e-20); + let p = f64::from_bits(1); + // pi * p would round to 3 * p here + approx::assert_relative_eq!(d.inverse_cdf(p), -6.442663821359281e302, max_relative = 1e-14); + let u = create_ok(0.0, 1.0); + approx::assert_relative_eq!(u.inverse_cdf(1e-12), -1.0 / (core::f64::consts::PI * 1e-12), max_relative = 1e-15); + } } From 01aab4723a0db057c590e3c1cdfe5eff1bdcade1 Mon Sep 17 00:00:00 2001 From: Dundy Pasupuleti Date: Sat, 31 Oct 2026 18:10:03 -0400 Subject: [PATCH 3/3] fix(cauchy): keep subnormal scale from underflowing in the lower tail Signed-off-by: Dundy Pasupuleti --- src/distribution/cauchy.rs | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/distribution/cauchy.rs b/src/distribution/cauchy.rs index ef433e2c..7360339a 100644 --- a/src/distribution/cauchy.rs +++ b/src/distribution/cauchy.rs @@ -171,10 +171,18 @@ impl ContinuousCDF for Cauchy { // every lower quantile below 2^-55 collapse to one value (gh #471) if x == 0.5 { self.location + } else if x == 0.0 { + f64::NEG_INFINITY } else if x < 1e-9 { // tan(pi x) = pi x to double precision here; dividing by x - // directly keeps a subnormal x from rounding pi x - self.location - self.scale * f64_consts::FRAC_1_PI / x + // directly keeps a subnormal x from rounding pi x; a + // subnormal scale is divided first so it does not underflow + let q = if self.scale < f64::MIN_POSITIVE { + self.scale / x * f64_consts::FRAC_1_PI + } else { + self.scale * f64_consts::FRAC_1_PI / x + }; + self.location - q } else if x < 0.5 { self.location - self.scale / (f64_consts::PI * x).tan() } else { @@ -555,6 +563,18 @@ mod tests { approx::assert_relative_eq!(d.inverse_cdf(1.0 - 1e-3), 318.30883898555044, max_relative = 1e-12); } + #[test] + fn test_inverse_cdf_subnormal_scale() { + let d = create_ok(0.0, f64::from_bits(1)); + approx::assert_relative_eq!( + d.inverse_cdf(1e-12), + -f64::from_bits(1) / (core::f64::consts::PI * 1e-12), + max_relative = 1e-3 + ); + assert!(d.inverse_cdf(1e-12) < 0.0); + assert_eq!(d.inverse_cdf(0.0), f64::NEG_INFINITY); + } + #[test] fn test_inverse_cdf_subnormal_p() { let d = create_ok(0.0, 1e-20);