diff --git a/src/macros/assert.rs b/src/macros/assert.rs
index 11c1d89d..b813a165 100644
--- a/src/macros/assert.rs
+++ b/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(
diff --git a/src/structs/dmat.rs b/src/structs/dmat.rs
index 757770a6..eec9cf5d 100644
--- a/src/structs/dmat.rs
+++ b/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> {
diff --git a/src/structs/dvec_macros.rs b/src/structs/dvec_macros.rs
index 8ee1f762..a8ef9bc9 100644
--- a/src/structs/dvec_macros.rs
+++ b/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> {
diff --git a/src/structs/iso_macros.rs b/src/structs/iso_macros.rs
index 98355a4b..4850e23e 100644
--- a/src/structs/iso_macros.rs
+++ b/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)
+            }
         }
     )
 );
diff --git a/src/structs/mat_macros.rs b/src/structs/mat_macros.rs
index 9d91d5e5..b9b0be5b 100644
--- a/src/structs/mat_macros.rs
+++ b/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))
+        }
     }
   )
 );
diff --git a/src/structs/quat.rs b/src/structs/quat.rs
index 4d4ef875..4f6472d5 100644
--- a/src/structs/quat.rs
+++ b/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> {
diff --git a/src/structs/rot_macros.rs b/src/structs/rot_macros.rs
index d6f6d81f..cfe331ec 100644
--- a/src/structs/rot_macros.rs
+++ b/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)
+            }
         }
     )
 );
diff --git a/src/structs/spec/vec0.rs b/src/structs/spec/vec0.rs
index 683b3e12..f44287a0 100644
--- a/src/structs/spec/vec0.rs
+++ b/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> {
diff --git a/src/structs/vec_macros.rs b/src/structs/vec_macros.rs
index f33e0644..44267b0c 100644
--- a/src/structs/vec_macros.rs
+++ b/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))*
+            }
         }
     )
 );
diff --git a/src/traits/operations.rs b/src/traits/operations.rs
index a1eac60c..64883e69 100644
--- a/src/traits/operations.rs
+++ b/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.