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

@@ -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::from(2.).unwrap();
+        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());
+    }
 }