diff --git a/src/macros/assert.rs b/src/macros/assert.rs
index 11c1d89d..de93983b 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 2ff5c362..4a58c9f6 100644
--- a/src/structs/dmat.rs
+++ b/src/structs/dmat.rs
@@ -605,11 +605,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 4c13b702..67eb5353 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 71d4c1e9..5100f48a 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 d71f71d3..8c926d5d 100644
--- a/src/structs/quat.rs
+++ b/src/structs/quat.rs
@@ -291,10 +291,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 2120ace4..8deaf8b0 100644
--- a/src/structs/rot_macros.rs
+++ b/src/structs/rot_macros.rs
@@ -250,6 +250,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)
@@ -259,6 +264,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 fe1607ed..bec5d088 100644
--- a/src/structs/spec/vec0.rs
+++ b/src/structs/spec/vec0.rs
@@ -205,10 +205,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 ae4dcad3..0d84234e 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 d731d13f..6f4461fd 100644
--- a/src/traits/operations.rs
+++ b/src/traits/operations.rs
@@ -158,6 +158,12 @@ pub trait ApproxEq<Eps>: Sized {
     /// 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 {
@@ -175,6 +181,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 {
@@ -187,6 +212,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.
diff --git a/tests/assert.rs b/tests/assert.rs
index 3b181aa6..9348a0f1 100644
--- a/tests/assert.rs
+++ b/tests/assert.rs
@@ -33,3 +33,36 @@ fn assert_approx_eq_eps_f32() {
 fn assert_approx_eq_eps_f64_fail() {
     assert_approx_eq_eps!(1.0f64, 1.1, 0.05);
 }
+
+#[test]
+fn assert_approx_eq_ulps_f32() {
+    let x = 1000000_f32;
+    let y = 1000000.1_f32;
+    assert!(x != y);
+    assert_approx_eq_ulps!(x, y, 3);
+}
+
+#[test]
+#[should_fail]
+fn assert_approx_eq_ulps_f32_fail() {
+    let x = 1000000_f32;
+    let y = 1000000.1_f32;
+    assert_approx_eq_ulps!(x, y, 2);
+}
+
+#[test]
+fn assert_approx_eq_ulps_f64() {
+    let x = 1000000_f64;
+    let y = 1000000.0000000003_f64;
+    assert!(x != y);
+    assert_approx_eq_ulps!(x, y, 4);
+}
+
+#[test]
+#[should_fail]
+fn assert_approx_eq_ulps_f64_fail() {
+    let x = 1000000_f64;
+    let y = 1000000.0000000003_f64;
+    assert!(x != y);
+    assert_approx_eq_ulps!(x, y, 3);
+}