From 86270c92ce0b1b27fd4c74ee46504c9136db6ffa Mon Sep 17 00:00:00 2001 From: Lukas Kalbertodt Date: Mon, 1 May 2023 18:35:32 +0200 Subject: [PATCH] Fix PERT distribution for when `mode` is very close to `(max - min) / 2` There is already a special condition for this case, as the general formula breaks for floats. However, the check uses `==` which is fishy for floats. This commit instead checks if the difference is smaller than the machine epsilon. Without this commit, this returns an error (despite being totally valid parameters for PERT): Pert::new(0.0, 0.48258883, 0.24129441) --- rand_distr/src/pert.rs | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/rand_distr/src/pert.rs b/rand_distr/src/pert.rs index db89fff7bfb..9ee44a3abb4 100644 --- a/rand_distr/src/pert.rs +++ b/rand_distr/src/pert.rs @@ -99,7 +99,14 @@ where let range = max - min; let mu = (min + max + shape * mode) / (shape + F::from(2.).unwrap()); - let v = if mu == mode { + + // We need to check if `mu == mode`, but due to rounding errors we can't + // use a direct comparison. The epsilon used for the comparison is + // twice the machine epsilon scaled by the larger of the two numbers. + // + // https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ + let epsilon = F::max(mu.abs(), mode.abs()) * F::epsilon() * (F::one() + F::one()); + let v = if (mu - mode).abs() < epsilon { shape * F::from(0.5).unwrap() + F::from(1.).unwrap() } else { (mu - min) * (F::from(2.).unwrap() * mode - min - max) / ((mode - mu) * (max - min)) @@ -151,4 +158,9 @@ mod test { fn pert_distributions_can_be_compared() { assert_eq!(Pert::new(1.0, 3.0, 2.0), Pert::new(1.0, 3.0, 2.0)); } + + #[test] + fn pert_mode_almost_half_range() { + assert!(Pert::new(0.0f32, 0.48258883, 0.24129441).is_ok()); + } }