ApproxEq trait enhanced with ULPs method of specifying closeness:

approx_eq_ulps() allows specification of epsilon as an integer number
  of Units in the Last Place (ULPs) difference between the two floating
  point values

  default approx_ulps() is set to 8.

  approx_eq() function continues to use epsilon method, although I
  recommend further commits and a migration towards the ULPs method.

src/macros/assert.rs

@@ -13,6 +13,21 @@ macro_rules! assert_approx_eq_eps(
     })
 );
 
+/// Asserts approximate equality within a given tolerance of two values with the
+/// `ApproxEq` trait, with tolerance specified in ULPs.
+#[macro_export]
+macro_rules! assert_approx_eq_ulps(
+    ($given: expr, $expected: expr, $ulps: expr) => ({
+        let ulps = &($ulps);
+        let (given_val, expected_val) = (&($given), &($expected));
+        if !ApproxEq::approx_eq_ulps(given_val, expected_val, ulps) {
+            panic!("assertion failed: `left ≈ right` (left: `{}`, right: `{}`, tolerance: `{}`)",
+                *given_val, *expected_val, *ulps
+            )
+        }
+    })
+);
+
 /// Asserts approximate equality of two values with the `ApproxEq` trait.
 #[macro_export]
 macro_rules! assert_approx_eq(

src/structs/dmat.rs

@@ -603,11 +603,22 @@ impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> {
         ApproxEq::approx_epsilon(None::<N>)
     }
 
+    #[inline]
+    fn approx_ulps(_: Option<DMat<N>>) -> u32 {
+        ApproxEq::approx_ulps(None::<N>)
+    }
+
     #[inline]
     fn approx_eq_eps(&self, other: &DMat<N>, epsilon: &N) -> bool {
         let zip = self.mij.iter().zip(other.mij.iter());
         zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon))
     }
+
+    #[inline]
+    fn approx_eq_ulps(&self, other: &DMat<N>, ulps: u32) -> bool {
+        let zip = self.mij.iter().zip(other.mij.iter());
+        zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps))
+    }
 }
 
 impl<N: Show + Copy> Show for DMat<N> {

src/structs/dvec_macros.rs

@@ -300,11 +300,22 @@ macro_rules! dvec_impl(
                 ApproxEq::approx_epsilon(None::<N>)
             }
 
+            #[inline]
+            fn approx_ulps(_: Option<$dvec<N>>) -> u32 {
+                ApproxEq::approx_ulps(None::<N>)
+            }
+
             #[inline]
             fn approx_eq_eps(&self, other: &$dvec<N>, epsilon: &N) -> bool {
                 let zip = self.as_slice().iter().zip(other.as_slice().iter());
                 zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon))
             }
+
+            #[inline]
+            fn approx_eq_ulps(&self, other: &$dvec<N>, ulps: u32) -> bool {
+                let zip = self.as_slice().iter().zip(other.as_slice().iter());
+                zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps))
+            }
         }
 
         impl<N: Copy + Mul<N, N> + Zero> Mul<N, $dvec<N>> for $dvec<N> {

src/structs/iso_macros.rs

@@ -316,11 +316,22 @@ macro_rules! approx_eq_impl(
                 ApproxEq::approx_epsilon(None::<N>)
             }
 
+            #[inline]
+            fn approx_ulps(_: Option<$t<N>>) -> u32 {
+                ApproxEq::approx_ulps(None::<N>)
+            }
+
             #[inline]
             fn approx_eq_eps(&self, other: &$t<N>, epsilon: &N) -> bool {
                 ApproxEq::approx_eq_eps(&self.rotation, &other.rotation, epsilon) &&
                     ApproxEq::approx_eq_eps(&self.translation, &other.translation, epsilon)
             }
+
+            #[inline]
+            fn approx_eq_ulps(&self, other: &$t<N>, ulps: u32) -> bool {
+                ApproxEq::approx_eq_ulps(&self.rotation, &other.rotation, ulps) &&
+                    ApproxEq::approx_eq_ulps(&self.translation, &other.translation, ulps)
+            }
         }
     )
 );

src/structs/mat_macros.rs

@@ -622,11 +622,22 @@ macro_rules! approx_eq_impl(
             ApproxEq::approx_epsilon(None::<N>)
         }
 
+        #[inline]
+        fn approx_ulps(_: Option<$t<N>>) -> u32 {
+            ApproxEq::approx_ulps(None::<N>)
+        }
+
         #[inline]
         fn approx_eq_eps(&self, other: &$t<N>, epsilon: &N) -> bool {
             let zip = self.iter().zip(other.iter());
             zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon))
         }
+
+        #[inline]
+        fn approx_eq_ulps(&self, other: &$t<N>, ulps: u32) -> bool {
+            let zip = self.iter().zip(other.iter());
+            zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps))
+        }
     }
   )
 );

src/structs/quat.rs

@@ -289,10 +289,20 @@ impl<N: ApproxEq<N>> ApproxEq<N> for UnitQuat<N> {
         ApproxEq::approx_epsilon(None::<N>)
     }
 
+    #[inline]
+    fn approx_ulps(_: Option<UnitQuat<N>>) -> u32 {
+        ApproxEq::approx_ulps(None::<N>)
+    }
+
     #[inline]
     fn approx_eq_eps(&self, other: &UnitQuat<N>, eps: &N) -> bool {
         ApproxEq::approx_eq_eps(&self.q, &other.q, eps)
     }
+
+    #[inline]
+    fn approx_eq_ulps(&self, other: &UnitQuat<N>, ulps: u32) -> bool {
+        ApproxEq::approx_eq_ulps(&self.q, &other.q, ulps)
+    }
 }
 
 impl<N: BaseFloat + ApproxEq<N>> Div<UnitQuat<N>, UnitQuat<N>> for UnitQuat<N> {

src/structs/rot_macros.rs

@@ -256,6 +256,11 @@ macro_rules! approx_eq_impl(
                 ApproxEq::approx_epsilon(None::<N>)
             }
 
+            #[inline]
+            fn approx_ulps(_: Option<$t<N>>) -> u32 {
+                ApproxEq::approx_ulps(None::<N>)
+            }
+
             #[inline]
             fn approx_eq(&self, other: &$t<N>) -> bool {
                 ApproxEq::approx_eq(&self.submat, &other.submat)
@@ -265,6 +270,11 @@ macro_rules! approx_eq_impl(
             fn approx_eq_eps(&self, other: &$t<N>, epsilon: &N) -> bool {
                 ApproxEq::approx_eq_eps(&self.submat, &other.submat, epsilon)
             }
+
+            #[inline]
+            fn approx_eq_ulps(&self, other: &$t<N>, ulps: u32) -> bool {
+                ApproxEq::approx_eq_ulps(&self.submat, &other.submat, ulps)
+            }
         }
     )
 );

src/structs/spec/vec0.rs

@@ -204,10 +204,19 @@ impl<N: ApproxEq<N>> ApproxEq<N> for vec::Vec0<N> {
         ApproxEq::approx_epsilon(None::<N>)
     }
 
+    fn approx_ulps(_: Option<vec::Vec0<N>>) -> u32 {
+        ApproxEq::approx_ulps(None::<N>)
+    }
+
     #[inline]
     fn approx_eq_eps(&self, _: &vec::Vec0<N>, _: &N) -> bool {
         true
     }
+
+    #[inline]
+    fn approx_eq_ulps(&self, _: &vec::Vec0<N>, _: u32) -> bool {
+        true
+    }
 }
 
 impl<N: One> One for vec::Vec0<N> {

src/structs/vec_macros.rs

@@ -609,6 +609,11 @@ macro_rules! approx_eq_impl(
                 ApproxEq::approx_epsilon(None::<N>)
             }
 
+            #[inline]
+            fn approx_ulps(_: Option<$t<N>>) -> u32 {
+                ApproxEq::approx_ulps(None::<N>)
+            }
+
             #[inline]
             fn approx_eq(&self, other: &$t<N>) -> bool {
                 ApproxEq::approx_eq(&self.$comp0, &other.$comp0)
@@ -620,6 +625,12 @@ macro_rules! approx_eq_impl(
                 ApproxEq::approx_eq_eps(&self.$comp0, &other.$comp0, eps)
                 $(&& ApproxEq::approx_eq_eps(&self.$compN, &other.$compN, eps))*
             }
+
+            #[inline]
+            fn approx_eq_ulps(&self, other: &$t<N>, ulps: u32) -> bool {
+                ApproxEq::approx_eq_ulps(&self.$comp0, &other.$comp0, ulps)
+                $(&& ApproxEq::approx_eq_ulps(&self.$compN, &other.$compN, ulps))*
+            }
         }
     )
 );

src/traits/operations.rs

@@ -156,6 +156,12 @@ pub trait ApproxEq<Eps> {
     /// Tests approximate equality using a custom epsilon.
     fn approx_eq_eps(&self, other: &Self, epsilon: &Eps) -> bool;
 
+    /// Default ULPs for approximation.
+    fn approx_ulps(unused_self: Option<Self>) -> u32;
+
+    /// Tests approximate equality using units in the last place (ULPs)
+    fn approx_eq_ulps(&self, other: &Self, ulps: u32) -> bool;
+
     /// Tests approximate equality.
     #[inline]
     fn approx_eq(&self, other: &Self) -> bool {
@@ -173,6 +179,25 @@ impl ApproxEq<f32> for f32 {
     fn approx_eq_eps(&self, other: &f32, epsilon: &f32) -> bool {
         ::abs(&(*self - *other)) < *epsilon
     }
+
+    fn approx_ulps(_: Option<f32>) -> u32 {
+        8
+    }
+
+    fn approx_eq_ulps(&self, other: &f32, ulps: u32) -> bool {
+        // Handle -0 == +0
+        if *self == *other { return true; }
+
+        // Otherwise, differing signs should be not-equal, even if within ulps
+        if self.signum() != other.signum() { return false; }
+
+        // IEEE754 floats are in the same order as 2s complement ints
+        // so this trick (subtracting the ints) works.
+        let iself: i32 = unsafe { ::std::mem::transmute(*self) };
+        let iother: i32 = unsafe { ::std::mem::transmute(*other) };
+
+        (iself - iother).abs() < ulps as i32
+    }
 }
 
 impl ApproxEq<f64> for f64 {
@@ -185,6 +210,23 @@ impl ApproxEq<f64> for f64 {
     fn approx_eq_eps(&self, other: &f64, approx_epsilon: &f64) -> bool {
         ::abs(&(*self - *other)) < *approx_epsilon
     }
+
+    fn approx_ulps(_: Option<f64>) -> u32 {
+        8
+    }
+
+    fn approx_eq_ulps(&self, other: &f64, ulps: u32) -> bool {
+        // Handle -0 == +0
+        if *self == *other { return true; }
+
+        // Otherwise, differing signs should be not-equal, even if within ulps
+        if self.signum() != other.signum() { return false; }
+
+        let iself: i64 = unsafe { ::std::mem::transmute(*self) };
+        let iother: i64 = unsafe { ::std::mem::transmute(*other) };
+
+        (iself - iother).abs() < ulps as i64
+    }
 }
 
 /// Trait of objects having an absolute value.