use another definition of ULP to compute ULP diff + extend explanation of this

Author: LorrensP-2158466

tests/pass/float.rs

@@ -12,15 +12,22 @@ use std::fmt::{Debug, Display, LowerHex};
 use std::hint::black_box;
 use std::{f32, f64};
 
-/// Another way of checking if 2 floating point numbers are almost equal to eachother.
-/// Using `a` and `b` as floating point numbers:
+/// Another way of checking if 2 floating-point numbers are almost equal to eachother.
+/// Using `a` and `b` as floating-point numbers:
 ///
-/// Instead of doing a simple EPSILON check (which we did at first):
-/// Absolute difference of `a` and `b` can't be greater than some E (10^-6, ...).
-/// `(a - b).asb() <= E`
+/// Instead of performing a simple EPSILON check (which we used at first):
+/// The absolute difference between 'a' and 'b' must not be greater than some E (10^-6, ...)
 ///
-/// We will now use ULP: `Units in the Last Place` or `Units of Least Precision`
-/// It is the absolute difference of the *unsigned integer* representation of `a` and `b`.
+/// We will now use ULP: `Units in the Last Place` or `Units of Least Precision`,
+/// more specific, the difference in ULP of `a` and `b`.
+/// First: The ULP of a float 'a' is the smallest possible change at 'a', so the ULP difference represents how
+/// many discrete floating-point steps are needed to reach 'b' from 'a'.
+///
+/// There are 2 techniques we can use to compute the ULP difference of 2 floating-point numbers:
+///
+/// ULP(a) is the difference between 'a' and the next larger representable floating-point number.
+/// So to get the ULP difference of `a` and `b`, we can take the absolute difference of
+/// the *unsigned integer* representation of `a` and `b`.
 /// For example checking 2 f32's, which in IEEE format looks like this:
 ///
 /// s: sign bit
@@ -29,71 +36,78 @@ use std::{f32, f64};
 /// f32: s | eeee eeee | mmm mmmm mmmm mmmm mmmm mmmm
 /// u32: seee eeee emmm mmmm mmmm mmmm mmmm mmmm
 ///
-/// We see that checking `a` and `b` with different signs has no meaning, but we should not forget
-/// -0.0 and +0.0.
-/// Same with exponents but no zero checking.
 ///
-/// So when Sign and Exponent are the same, we can get a reasonable ULP value
-/// by doing the absolute difference, and if this is less than or equal to our chosen upper bound
-/// we can say that `a` and `b` are approximately equal.
+/// So when the sign and exponent are the same, we can compute a reasonable ULP difference.
+///
+/// But if you know some details of floating-point numbers, it is possible to get 2 values that sit on the border of the exponent,
+/// so `a` can have exponent but `b`  can have a different exponent while still being approximately equal to `a`. For this we can use
+/// John Harrison's definition. Which states that ULP(a) is the distance between the 2 closest floating-point numbers
+/// `x` and `y` around `a`, satisfying x < a < y, x != y. This makes the ULP of `a` twice as large as in the previous defintion and so
+/// if we want to use this, we have to halve it.
+/// This ULP value of `a` can then be used to find the ULP difference of `a` and `b`:
+/// Take the difference of `b` and `a` and divide it by that ULP and finally round it. This is the same value as the first definition,
+/// but this way can go around that exponent border problem.
+///
+/// If this ULP difference is less than or equal to our chosen upper bound
+/// we can say that `a` and `b` are approximately equal, because they lie "close" enough to each other to be considered equal.
 ///
-/// Basically ULP can be seen as a distance metric of floating point numbers, but having
-/// the same amount of "spacing" between all consecutive representable values. So eventhough 2 very large floating point numbers
+/// Note: We can see that checking `a` and `b` with different signs has no meaning, but we should not forget
+/// -0.0 and +0.0.
+///
+/// Essentially ULP can be seen as a distance metric of floating-point numbers, but with
+/// the same amount of "spacing" between all consecutive representable values. So even though 2 very large floating point numbers
 /// have a large value difference, their ULP can still be 1, so they are still "approximatly equal",
 /// but the EPSILON check would have failed.
 fn approx_eq_check<F: Float>(
     actual: F,
     expected: F,
     allowed_ulp: F::Int,
-) -> Result<(), NotApproxEq<F::Int>>
+) -> Result<(), NotApproxEq<F>>
 where
-    F::Int: PartialOrd + AbsDiff,
+    F::Int: PartialOrd,
 {
     let actual_signum = actual.signum();
     let expected_signum = expected.signum();
 
     if actual_signum != expected_signum {
+        // Floats with different signs must both be 0.
         if actual != expected {
             return Err(NotApproxEq::SignsDiffer);
         }
     } else {
-        // reinterpret the floating point numbers as unsigned integers
-        let actual_bits = actual.to_bits();
-        let expected_bits = expected.to_bits();
+        let ulp = (expected.next_up() - expected.next_down()).halve();
+        let ulp_diff = ((actual - expected) / ulp).round().as_int();
 
-        // take their absolute difference to get ULP
-        let ulp_diff = actual_bits.abs_diff(expected_bits);
         if ulp_diff > allowed_ulp {
             return Err(NotApproxEq::UlpFail(ulp_diff));
         }
     }
     Ok(())
 }
 
-/// Give more context to execution and result of [`approx_eq_check`]
-enum NotApproxEq<Int> {
+/// Give more context to execution and result of [`approx_eq_check`].
+enum NotApproxEq<F: Float> {
     SignsDiffer,
 
-    /// Contains the actual ulp value calculated
-    UlpFail(Int),
+    /// Contains the actual ulp value calculated.
+    UlpFail(F::Int),
 }
 
 macro_rules! assert_approx_eq {
     ($a:expr, $b:expr, $ulp:expr) => {{
         let (a, b) = ($a, $b);
-        let ulp = $ulp;
-        // Floats with different signs must both be 0.
-        match approx_eq_check(a, b, ulp) {
+        let allowed_ulp = $ulp;
+        match approx_eq_check(a, b, allowed_ulp) {
             Err(NotApproxEq::SignsDiffer) =>
                 panic!("{a:?} is not approximately equal to {b:?}: signs differ"),
             Err(NotApproxEq::UlpFail(actual_ulp)) =>
-                panic!("{a:?} is not approximately equal to {b:?}\nulp diff: {actual_ulp} > {ulp}"),
-            _ => {}
+                panic!("{a:?} is not approximately equal to {b:?}\nulp diff: {actual_ulp} > {allowed_ulp}"),
+            Ok(_) => {}
         };
     }};
 
     ($a:expr, $b: expr) => {
-        // accept up to 64ULP (16ULP for host floats and 16ULP for artificial error and 32 for any rounding errors)
+        // accept up to 64ULP (16ULP for host floats and 16ULP for miri artificial error and 32 for any rounding errors)
         assert_approx_eq!($a, $b, 64);
     };
 }
@@ -114,12 +128,14 @@ fn main() {
     test_non_determinism();
 }
 
-// to help `approx_eq_check`
-trait AbsDiff {
-    fn abs_diff(self, other: Self) -> Self;
-}
-
-trait Float: Copy + PartialEq + Debug {
+trait Float:
+    Copy
+    + PartialEq
+    + Debug
+    + std::ops::Sub<Output = Self>
+    + std::cmp::PartialOrd
+    + std::ops::Div<Output = Self>
+{
     /// The unsigned integer with the same bit width as this float
     type Int: Copy + PartialEq + LowerHex + Debug;
     const BITS: u32 = size_of::<Self>() as u32 * 8;
@@ -134,7 +150,14 @@ trait Float: Copy + PartialEq + Debug {
 
     fn to_bits(self) -> Self::Int;
 
+    // to make "approx_eq_check" generic
     fn signum(self) -> Self;
+    fn next_up(self) -> Self;
+    fn next_down(self) -> Self;
+    fn round(self) -> Self;
+    // self / 2
+    fn halve(self) -> Self;
+    fn as_int(self) -> Self::Int;
 }
 
 macro_rules! impl_float {
@@ -151,11 +174,22 @@ macro_rules! impl_float {
             fn signum(self) -> Self {
                 self.signum()
             }
-        }
+            fn next_up(self) -> Self {
+                self.next_up()
+            }
+            fn next_down(self) -> Self {
+                self.next_down()
+            }
+            fn round(self) -> Self {
+                self.round()
+            }
+
+            fn halve(self) -> Self {
+                self / 2.0
+            }
 
-        impl AbsDiff for $ity {
-            fn abs_diff(self, other: Self) -> Self {
-                self.abs_diff(other)
+            fn as_int(self) -> Self::Int {
+                self as Self::Int
             }
         }
     };