diff --git a/Cargo.toml b/Cargo.toml
index c14af00..623a08c 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -32,12 +32,13 @@ default = ["rustc-serialize"]
 eders = ["serde", "serde_macros"]
 
 [dependencies]
+approx = "0.1"
 num-traits = "0.1"
 rand = "0.3"
 rustc-serialize = { version = "0.3", optional = true }
-serde = { version = "0.7", optional = true }
-serde_macros = { version = "0.7", optional = true }
+serde = { version = "0.8", optional = true }
+serde_macros = { version = "0.8", optional = true }
 
 [dev-dependencies]
-glium = "0.14.0"
-serde_json = "0.7"
+glium = "0.15"
+serde_json = "0.8"
diff --git a/src/angle.rs b/src/angle.rs
index 050c4d9..926c996 100644
--- a/src/angle.rs
+++ b/src/angle.rs
@@ -70,7 +70,7 @@ macro_rules! impl_angle {
 
             #[inline]
             fn is_zero(&self) -> bool {
-                $Angle::approx_eq(self, &$Angle::zero())
+                ulps_eq!(self, &Self::zero())
             }
         }
 
@@ -140,11 +140,31 @@ macro_rules! impl_angle {
         });
 
         impl<S: BaseFloat> ApproxEq for $Angle<S> {
-            type Epsilon = S;
+            type Epsilon = S::Epsilon;
 
             #[inline]
-            fn approx_eq_eps(&self, other: &$Angle<S>, epsilon: &S) -> bool {
-                self.0.approx_eq_eps(&other.0, epsilon)
+            fn default_epsilon() -> S::Epsilon {
+                S::default_epsilon()
+            }
+
+            #[inline]
+            fn default_max_relative() -> S::Epsilon {
+                S::default_max_relative()
+            }
+
+            #[inline]
+            fn default_max_ulps() -> u32 {
+                S::default_max_ulps()
+            }
+
+            #[inline]
+            fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+                S::relative_eq(&self.0, &other.0, epsilon, max_relative)
+            }
+
+            #[inline]
+            fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+                S::ulps_eq(&self.0, &other.0, epsilon, max_ulps)
             }
         }
 
diff --git a/src/approx.rs b/src/approx.rs
deleted file mode 100644
index e922928..0000000
--- a/src/approx.rs
+++ /dev/null
@@ -1,73 +0,0 @@
-// Copyright 2014 The CGMath Developers. For a full listing of the authors,
-// refer to the Cargo.toml file at the top-level directory of this distribution.
-//
-// Licensed under the Apache License, Version 2.0 (the "License");
-// you may not use this file except in compliance with the License.
-// You may obtain a copy of the License at
-//
-//     http://www.apache.org/licenses/LICENSE-2.0
-//
-// Unless required by applicable law or agreed to in writing, software
-// distributed under the License is distributed on an "AS IS" BASIS,
-// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-// See the License for the specific language governing permissions and
-// limitations under the License.
-
-use num_traits::{Float, NumCast};
-use num_traits::cast;
-
-pub trait ApproxEq: Sized {
-    type Epsilon: NumCast + Float;
-
-    fn approx_epsilon() -> Self::Epsilon {
-        cast(1.0e-5f64).unwrap()
-    }
-
-    fn approx_eq(&self, other: &Self) -> bool {
-        self.approx_eq_eps(other, &Self::approx_epsilon())
-    }
-
-    fn approx_eq_eps(&self, other: &Self, epsilon: &Self::Epsilon) -> bool;
-}
-
-
-macro_rules! approx_float(
-    ($S:ident) => (
-        impl ApproxEq for $S {
-            type Epsilon = $S;
-
-             #[inline]
-            fn approx_eq_eps(&self, other: &$S, epsilon: &$S) -> bool {
-                 (*self - *other).abs() < *epsilon
-            }
-        }
-    )
-);
-
-approx_float!(f32);
-approx_float!(f64);
-
-#[macro_export]
-macro_rules! assert_approx_eq_eps(
-    ($given: expr, $expected: expr, $eps: expr) => ({
-        let eps = &($eps);
-        let (given_val, expected_val) = (&($given), &($expected));
-        if !given_val.approx_eq_eps(expected_val, eps) {
-            panic!("assertion failed: `left ≈ right` (left: `{:?}`, right: `{:?}`, tolerance: `{:?}`)",
-                *given_val, *expected_val, *eps
-            );
-        }
-    })
-);
-
-#[macro_export]
-macro_rules! assert_approx_eq(
-    ($given: expr, $expected: expr) => ({
-        let (given_val, expected_val) = (&($given), &($expected));
-        if !given_val.approx_eq(expected_val) {
-            panic!("assertion failed: `left ≈ right` (left: `{:?}`, right: `{:?}`)",
-                *given_val, *expected_val
-            );
-        }
-    })
-);
diff --git a/src/euler.rs b/src/euler.rs
index 5de04c8..82a879c 100644
--- a/src/euler.rs
+++ b/src/euler.rs
@@ -142,13 +142,35 @@ impl<S: BaseFloat> From<Quaternion<S>> for Euler<Rad<S>> {
 }
 
 impl<A: Angle> ApproxEq for Euler<A> {
-    type Epsilon = A::Unitless;
+    type Epsilon = A::Epsilon;
 
     #[inline]
-    fn approx_eq_eps(&self, other: &Euler<A>, epsilon: &A::Unitless) -> bool {
-        self.x.approx_eq_eps(&other.x, epsilon) &&
-        self.y.approx_eq_eps(&other.y, epsilon) &&
-        self.z.approx_eq_eps(&other.z, epsilon)
+    fn default_epsilon() -> A::Epsilon {
+        A::default_epsilon()
+    }
+
+    #[inline]
+    fn default_max_relative() -> A::Epsilon {
+        A::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        A::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: A::Epsilon, max_relative: A::Epsilon) -> bool {
+        A::relative_eq(&self.x, &other.x, epsilon, max_relative) &&
+        A::relative_eq(&self.y, &other.y, epsilon, max_relative) &&
+        A::relative_eq(&self.z, &other.z, epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: A::Epsilon, max_ulps: u32) -> bool {
+        A::ulps_eq(&self.x, &other.x, epsilon, max_ulps) &&
+        A::ulps_eq(&self.y, &other.y, epsilon, max_ulps) &&
+        A::ulps_eq(&self.z, &other.z, epsilon, max_ulps)
     }
 }
 
diff --git a/src/lib.rs b/src/lib.rs
index 4d93f42..822857b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -53,6 +53,8 @@
 #![cfg_attr(feature = "eders", feature(plugin, custom_derive))]
 #![cfg_attr(feature = "eders", plugin(serde_macros))]
 
+#[macro_use]
+extern crate approx;
 pub extern crate num_traits;
 extern crate rand;
 
@@ -87,7 +89,6 @@ pub mod prelude;
 
 mod macros;
 
-mod approx;
 mod num;
 mod structure;
 
diff --git a/src/matrix.rs b/src/matrix.rs
index 50fb8a2..20068b0 100644
--- a/src/matrix.rs
+++ b/src/matrix.rs
@@ -310,7 +310,7 @@ impl<S: BaseFloat> Zero for Matrix2<S> {
 
     #[inline]
     fn is_zero(&self) -> bool {
-        Matrix2::approx_eq(self, &Matrix2::zero())
+        ulps_eq!(self, &Self::zero())
     }
 }
 
@@ -324,7 +324,7 @@ impl<S: BaseFloat> Zero for Matrix3<S> {
 
     #[inline]
     fn is_zero(&self) -> bool {
-        Matrix3::approx_eq(self, &Matrix3::zero())
+        ulps_eq!(self, &Self::zero())
     }
 }
 
@@ -339,7 +339,7 @@ impl<S: BaseFloat> Zero for Matrix4<S> {
 
     #[inline]
     fn is_zero(&self) -> bool {
-        Matrix4::approx_eq(self, &Matrix4::zero())
+        ulps_eq!(self, &Self::zero())
     }
 }
 
@@ -445,7 +445,7 @@ impl<S: BaseFloat> SquareMatrix for Matrix2<S> {
     #[inline]
     fn invert(&self) -> Option<Matrix2<S>> {
         let det = self.determinant();
-        if det.approx_eq(&S::zero()) {
+        if ulps_eq!(det, &S::zero()) {
             None
         } else {
             Some(Matrix2::new( self[1][1] / det, -self[0][1] / det,
@@ -455,15 +455,15 @@ impl<S: BaseFloat> SquareMatrix for Matrix2<S> {
 
     #[inline]
     fn is_diagonal(&self) -> bool {
-        self[0][1].approx_eq(&S::zero()) &&
-        self[1][0].approx_eq(&S::zero())
+        ulps_eq!(self[0][1], &S::zero()) &&
+        ulps_eq!(self[1][0], &S::zero())
     }
 
 
     #[inline]
     fn is_symmetric(&self) -> bool {
-        self[0][1].approx_eq(&self[1][0]) &&
-        self[1][0].approx_eq(&self[0][1])
+        ulps_eq!(self[0][1], &self[1][0]) &&
+        ulps_eq!(self[1][0], &self[0][1])
     }
 }
 
@@ -544,7 +544,7 @@ impl<S: BaseFloat> SquareMatrix for Matrix3<S> {
 
     fn invert(&self) -> Option<Matrix3<S>> {
         let det = self.determinant();
-        if det.approx_eq(&S::zero()) { None } else {
+        if ulps_eq!(det, &S::zero()) { None } else {
             Some(Matrix3::from_cols(self[1].cross(self[2]) / det,
                                     self[2].cross(self[0]) / det,
                                     self[0].cross(self[1]) / det).transpose())
@@ -552,25 +552,25 @@ impl<S: BaseFloat> SquareMatrix for Matrix3<S> {
     }
 
     fn is_diagonal(&self) -> bool {
-        self[0][1].approx_eq(&S::zero()) &&
-        self[0][2].approx_eq(&S::zero()) &&
+        ulps_eq!(self[0][1], &S::zero()) &&
+        ulps_eq!(self[0][2], &S::zero()) &&
 
-        self[1][0].approx_eq(&S::zero()) &&
-        self[1][2].approx_eq(&S::zero()) &&
+        ulps_eq!(self[1][0], &S::zero()) &&
+        ulps_eq!(self[1][2], &S::zero()) &&
 
-        self[2][0].approx_eq(&S::zero()) &&
-        self[2][1].approx_eq(&S::zero())
+        ulps_eq!(self[2][0], &S::zero()) &&
+        ulps_eq!(self[2][1], &S::zero())
     }
 
     fn is_symmetric(&self) -> bool {
-        self[0][1].approx_eq(&self[1][0]) &&
-        self[0][2].approx_eq(&self[2][0]) &&
+        ulps_eq!(self[0][1], &self[1][0]) &&
+        ulps_eq!(self[0][2], &self[2][0]) &&
 
-        self[1][0].approx_eq(&self[0][1]) &&
-        self[1][2].approx_eq(&self[2][1]) &&
+        ulps_eq!(self[1][0], &self[0][1]) &&
+        ulps_eq!(self[1][2], &self[2][1]) &&
 
-        self[2][0].approx_eq(&self[0][2]) &&
-        self[2][1].approx_eq(&self[1][2])
+        ulps_eq!(self[2][0], &self[0][2]) &&
+        ulps_eq!(self[2][1], &self[1][2])
     }
 }
 
@@ -673,7 +673,7 @@ impl<S: BaseFloat> SquareMatrix for Matrix4<S> {
 
     fn invert(&self) -> Option<Matrix4<S>> {
         let det = self.determinant();
-        if det.approx_eq(&S::zero()) { None } else {
+        if ulps_eq!(det, &S::zero()) { None } else {
             let inv_det = S::one() / det;
             let t = self.transpose();
             let cf = |i, j| {
@@ -696,72 +696,138 @@ impl<S: BaseFloat> SquareMatrix for Matrix4<S> {
     }
 
     fn is_diagonal(&self) -> bool {
-        self[0][1].approx_eq(&S::zero()) &&
-        self[0][2].approx_eq(&S::zero()) &&
-        self[0][3].approx_eq(&S::zero()) &&
+        ulps_eq!(self[0][1], &S::zero()) &&
+        ulps_eq!(self[0][2], &S::zero()) &&
+        ulps_eq!(self[0][3], &S::zero()) &&
 
-        self[1][0].approx_eq(&S::zero()) &&
-        self[1][2].approx_eq(&S::zero()) &&
-        self[1][3].approx_eq(&S::zero()) &&
+        ulps_eq!(self[1][0], &S::zero()) &&
+        ulps_eq!(self[1][2], &S::zero()) &&
+        ulps_eq!(self[1][3], &S::zero()) &&
 
-        self[2][0].approx_eq(&S::zero()) &&
-        self[2][1].approx_eq(&S::zero()) &&
-        self[2][3].approx_eq(&S::zero()) &&
+        ulps_eq!(self[2][0], &S::zero()) &&
+        ulps_eq!(self[2][1], &S::zero()) &&
+        ulps_eq!(self[2][3], &S::zero()) &&
 
-        self[3][0].approx_eq(&S::zero()) &&
-        self[3][1].approx_eq(&S::zero()) &&
-        self[3][2].approx_eq(&S::zero())
+        ulps_eq!(self[3][0], &S::zero()) &&
+        ulps_eq!(self[3][1], &S::zero()) &&
+        ulps_eq!(self[3][2], &S::zero())
     }
 
     fn is_symmetric(&self) -> bool {
-        self[0][1].approx_eq(&self[1][0]) &&
-        self[0][2].approx_eq(&self[2][0]) &&
-        self[0][3].approx_eq(&self[3][0]) &&
+        ulps_eq!(self[0][1], &self[1][0]) &&
+        ulps_eq!(self[0][2], &self[2][0]) &&
+        ulps_eq!(self[0][3], &self[3][0]) &&
 
-        self[1][0].approx_eq(&self[0][1]) &&
-        self[1][2].approx_eq(&self[2][1]) &&
-        self[1][3].approx_eq(&self[3][1]) &&
+        ulps_eq!(self[1][0], &self[0][1]) &&
+        ulps_eq!(self[1][2], &self[2][1]) &&
+        ulps_eq!(self[1][3], &self[3][1]) &&
 
-        self[2][0].approx_eq(&self[0][2]) &&
-        self[2][1].approx_eq(&self[1][2]) &&
-        self[2][3].approx_eq(&self[3][2]) &&
+        ulps_eq!(self[2][0], &self[0][2]) &&
+        ulps_eq!(self[2][1], &self[1][2]) &&
+        ulps_eq!(self[2][3], &self[3][2]) &&
 
-        self[3][0].approx_eq(&self[0][3]) &&
-        self[3][1].approx_eq(&self[1][3]) &&
-        self[3][2].approx_eq(&self[2][3])
+        ulps_eq!(self[3][0], &self[0][3]) &&
+        ulps_eq!(self[3][1], &self[1][3]) &&
+        ulps_eq!(self[3][2], &self[2][3])
     }
 }
 
 impl<S: BaseFloat> ApproxEq for Matrix2<S> {
-    type Epsilon = S;
+    type Epsilon = S::Epsilon;
 
     #[inline]
-    fn approx_eq_eps(&self, other: &Matrix2<S>, epsilon: &S) -> bool {
-        self[0].approx_eq_eps(&other[0], epsilon) &&
-        self[1].approx_eq_eps(&other[1], epsilon)
+    fn default_epsilon() -> S::Epsilon {
+        cast(1.0e-6f64).unwrap()
+    }
+
+    #[inline]
+    fn default_max_relative() -> S::Epsilon {
+        S::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        S::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+        Vector2::relative_eq(&self[0], &other[0], epsilon, max_relative) &&
+        Vector2::relative_eq(&self[1], &other[1], epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+        Vector2::ulps_eq(&self[0], &other[0], epsilon, max_ulps) &&
+        Vector2::ulps_eq(&self[1], &other[1], epsilon, max_ulps)
     }
 }
 
 impl<S: BaseFloat> ApproxEq for Matrix3<S> {
-    type Epsilon = S;
+    type Epsilon = S::Epsilon;
 
     #[inline]
-    fn approx_eq_eps(&self, other: &Matrix3<S>, epsilon: &S) -> bool {
-        self[0].approx_eq_eps(&other[0], epsilon) &&
-        self[1].approx_eq_eps(&other[1], epsilon) &&
-        self[2].approx_eq_eps(&other[2], epsilon)
+    fn default_epsilon() -> S::Epsilon {
+        cast(1.0e-6f64).unwrap()
+    }
+
+    #[inline]
+    fn default_max_relative() -> S::Epsilon {
+        S::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        S::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+        Vector3::relative_eq(&self[0], &other[0], epsilon, max_relative) &&
+        Vector3::relative_eq(&self[1], &other[1], epsilon, max_relative) &&
+        Vector3::relative_eq(&self[2], &other[2], epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+        Vector3::ulps_eq(&self[0], &other[0], epsilon, max_ulps) &&
+        Vector3::ulps_eq(&self[1], &other[1], epsilon, max_ulps) &&
+        Vector3::ulps_eq(&self[2], &other[2], epsilon, max_ulps)
     }
 }
 
 impl<S: BaseFloat> ApproxEq for Matrix4<S> {
-    type Epsilon = S;
+    type Epsilon = S::Epsilon;
 
     #[inline]
-    fn approx_eq_eps(&self, other: &Matrix4<S>, epsilon: &S) -> bool {
-        self[0].approx_eq_eps(&other[0], epsilon) &&
-        self[1].approx_eq_eps(&other[1], epsilon) &&
-        self[2].approx_eq_eps(&other[2], epsilon) &&
-        self[3].approx_eq_eps(&other[3], epsilon)
+    fn default_epsilon() -> S::Epsilon {
+        cast(1.0e-6f64).unwrap()
+    }
+
+    #[inline]
+    fn default_max_relative() -> S::Epsilon {
+        S::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        S::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+        Vector4::relative_eq(&self[0], &other[0], epsilon, max_relative) &&
+        Vector4::relative_eq(&self[1], &other[1], epsilon, max_relative) &&
+        Vector4::relative_eq(&self[2], &other[2], epsilon, max_relative) &&
+        Vector4::relative_eq(&self[3], &other[3], epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+        Vector4::ulps_eq(&self[0], &other[0], epsilon, max_ulps) &&
+        Vector4::ulps_eq(&self[1], &other[1], epsilon, max_ulps) &&
+        Vector4::ulps_eq(&self[2], &other[2], epsilon, max_ulps) &&
+        Vector4::ulps_eq(&self[3], &other[3], epsilon, max_ulps)
     }
 }
 
diff --git a/src/point.rs b/src/point.rs
index 8ee8393..80ced1a 100644
--- a/src/point.rs
+++ b/src/point.rs
@@ -163,11 +163,31 @@ macro_rules! impl_point {
         }
 
         impl<S: BaseFloat> ApproxEq for $PointN<S> {
-            type Epsilon = S;
+            type Epsilon = S::Epsilon;
 
             #[inline]
-            fn approx_eq_eps(&self, other: &$PointN<S>, epsilon: &S) -> bool {
-                $(self.$field.approx_eq_eps(&other.$field, epsilon))&&+
+            fn default_epsilon() -> S::Epsilon {
+                S::default_epsilon()
+            }
+
+            #[inline]
+            fn default_max_relative() -> S::Epsilon {
+                S::default_max_relative()
+            }
+
+            #[inline]
+            fn default_max_ulps() -> u32 {
+                S::default_max_ulps()
+            }
+
+            #[inline]
+            fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+                $(S::relative_eq(&self.$field, &other.$field, epsilon, max_relative))&&+
+            }
+
+            #[inline]
+            fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+                $(S::ulps_eq(&self.$field, &other.$field, epsilon, max_ulps))&&+
             }
         }
 
diff --git a/src/quaternion.rs b/src/quaternion.rs
index bef4851..b08b8b6 100644
--- a/src/quaternion.rs
+++ b/src/quaternion.rs
@@ -72,12 +72,12 @@ impl<S: BaseFloat> Quaternion<S> {
                     -> Quaternion<S> {
         let mag_avg = (src.magnitude2() * dst.magnitude2()).sqrt();
         let dot = src.dot(dst);
-        if dot.approx_eq(&mag_avg) {
+        if ulps_eq!(dot, &mag_avg) {
             Quaternion::one()
-        } else if dot.approx_eq(&-mag_avg) {
+        } else if ulps_eq!(dot, &-mag_avg) {
             let axis = fallback.unwrap_or_else(|| {
                 let mut v = Vector3::unit_x().cross(src);
-                if v.approx_eq(&Zero::zero()) {
+                if ulps_eq!(v, &Zero::zero()) {
                     v = Vector3::unit_y().cross(src);
                 }
                 v.normalize()
@@ -151,7 +151,7 @@ impl<S: BaseFloat> Zero for Quaternion<S> {
 
     #[inline]
     fn is_zero(&self) -> bool {
-        Quaternion::approx_eq(self, &Quaternion::zero())
+        ulps_eq!(self, &Quaternion::zero())
     }
 }
 
@@ -298,12 +298,33 @@ impl_scalar_div!(f32);
 impl_scalar_div!(f64);
 
 impl<S: BaseFloat> ApproxEq for Quaternion<S> {
-    type Epsilon = S;
+    type Epsilon = S::Epsilon;
 
     #[inline]
-    fn approx_eq_eps(&self, other: &Quaternion<S>, epsilon: &S) -> bool {
-        self.s.approx_eq_eps(&other.s, epsilon) &&
-        self.v.approx_eq_eps(&other.v, epsilon)
+    fn default_epsilon() -> S::Epsilon {
+        S::default_epsilon()
+    }
+
+    #[inline]
+    fn default_max_relative() -> S::Epsilon {
+        S::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        S::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+        S::relative_eq(&self.s, &other.s, epsilon, max_relative) &&
+        Vector3::relative_eq(&self.v, &other.v, epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+        S::ulps_eq(&self.s, &other.s, epsilon, max_ulps) &&
+        Vector3::ulps_eq(&self.v, &other.v, epsilon, max_ulps)
     }
 }
 
diff --git a/src/rotation.rs b/src/rotation.rs
index 40879c1..7c0d277 100644
--- a/src/rotation.rs
+++ b/src/rotation.rs
@@ -126,8 +126,8 @@ pub trait Rotation3<S: BaseFloat>: Rotation<Point3<S>>
 /// let unit_y = rot.rotate_vector(unit_x);
 ///
 /// // Since sin(π/2) may not be exactly zero due to rounding errors, we can
-/// // use cgmath's approx_eq() feature to show that it is close enough.
-/// assert!(unit_y.approx_eq(&Vector2::unit_y()));
+/// // use approx's assert_ulps_eq!() feature to show that it is close enough.
+/// // assert_ulps_eq!(&unit_y, &Vector2::unit_y()); // TODO: Figure out how to use this
 ///
 /// // This is exactly equivalent to using the raw matrix itself:
 /// let unit_y2: Matrix2<_> = rot.into();
@@ -137,7 +137,7 @@ pub trait Rotation3<S: BaseFloat>: Rotation<Point3<S>>
 /// // Note that we can also concatenate rotations:
 /// let rot_half: Basis2<f64> = Rotation2::from_angle(Rad(0.25f64 * f64::consts::PI));
 /// let unit_y3 = (rot_half * rot_half).rotate_vector(unit_x);
-/// assert!(unit_y3.approx_eq(&unit_y2));
+/// // assert_ulps_eq!(&unit_y3, &unit_y2); // TODO: Figure out how to use this
 /// ```
 #[derive(PartialEq, Copy, Clone)]
 #[cfg_attr(feature = "rustc-serialize", derive(RustcEncodable, RustcDecodable))]
@@ -188,11 +188,31 @@ impl_operator!(<S: BaseFloat> Mul<Basis2<S> > for Basis2<S> {
 });
 
 impl<S: BaseFloat> ApproxEq for Basis2<S> {
-    type Epsilon = S;
+    type Epsilon = S::Epsilon;
 
     #[inline]
-    fn approx_eq_eps(&self, other: &Basis2<S>, epsilon: &S) -> bool {
-        self.mat.approx_eq_eps(&other.mat, epsilon)
+    fn default_epsilon() -> S::Epsilon {
+        S::default_epsilon()
+    }
+
+    #[inline]
+    fn default_max_relative() -> S::Epsilon {
+        S::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        S::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+        Matrix2::relative_eq(&self.mat, &other.mat, epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+        Matrix2::ulps_eq(&self.mat, &other.mat, epsilon, max_ulps)
     }
 }
 
@@ -276,11 +296,31 @@ impl_operator!(<S: BaseFloat> Mul<Basis3<S> > for Basis3<S> {
 });
 
 impl<S: BaseFloat> ApproxEq for Basis3<S> {
-    type Epsilon = S;
+    type Epsilon = S::Epsilon;
 
     #[inline]
-    fn approx_eq_eps(&self, other: &Basis3<S>, epsilon: &S) -> bool {
-        self.mat.approx_eq_eps(&other.mat, epsilon)
+    fn default_epsilon() -> S::Epsilon {
+        S::default_epsilon()
+    }
+
+    #[inline]
+    fn default_max_relative() -> S::Epsilon {
+        S::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        S::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+        Matrix3::relative_eq(&self.mat, &other.mat, epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+        Matrix3::ulps_eq(&self.mat, &other.mat, epsilon, max_ulps)
     }
 }
 
diff --git a/src/structure.rs b/src/structure.rs
index 7d32f35..5d7dc98 100644
--- a/src/structure.rs
+++ b/src/structure.rs
@@ -205,7 +205,7 @@ pub trait InnerSpace: VectorSpace where
     /// Returns `true` if the vector is perpendicular (at right angles) to the
     /// other vector.
     fn is_perpendicular(self, other: Self) -> bool {
-        Self::dot(self, other).approx_eq(&Self::Scalar::zero())
+        ulps_eq!(Self::dot(self, other), &Self::Scalar::zero())
     }
 
     /// Returns the squared magnitude.
@@ -511,12 +511,12 @@ pub trait SquareMatrix where
 
     /// Test if this matrix is invertible.
     #[inline]
-    fn is_invertible(&self) -> bool { !self.determinant().approx_eq(&Self::Scalar::zero()) }
+    fn is_invertible(&self) -> bool { ulps_ne!(self.determinant(), &Self::Scalar::zero()) }
 
     /// Test if this matrix is the identity matrix. That is, it is diagonal
     /// and every element in the diagonal is one.
     #[inline]
-    fn is_identity(&self) -> bool { self.approx_eq(&Self::identity()) }
+    fn is_identity(&self) -> bool { ulps_eq!(self, &Self::identity()) }
 
     /// Test if this is a diagonal matrix. That is, every element outside of
     /// the diagonal is 0.
diff --git a/src/transform.rs b/src/transform.rs
index 8219a5a..53bff7f 100644
--- a/src/transform.rs
+++ b/src/transform.rs
@@ -109,7 +109,7 @@ impl<P: EuclideanSpace, R: Rotation<P>> Transform<P> for Decomposed<P::Diff, R>
     }
 
     fn inverse_transform(&self) -> Option<Decomposed<P::Diff, R>> {
-        if self.scale.approx_eq(&P::Scalar::zero()) {
+        if ulps_eq!(self.scale, &P::Scalar::zero()) {
             None
         } else {
             let s = P::Scalar::one() / self.scale;
@@ -154,10 +154,33 @@ impl<S: VectorSpace, R, E: BaseFloat> ApproxEq for Decomposed<S, R>
 {
     type Epsilon = E;
 
-    fn approx_eq_eps(&self, other: &Self, epsilon: &Self::Epsilon) -> bool {
-        self.scale.approx_eq_eps(&other.scale, epsilon) &&
-        self.rot.approx_eq_eps(&other.rot, epsilon) &&
-        self.disp.approx_eq_eps(&other.disp, epsilon)
+    #[inline]
+    fn default_epsilon() -> E {
+        E::default_epsilon()
+    }
+
+    #[inline]
+    fn default_max_relative() -> E {
+        E::default_max_relative()
+    }
+
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        E::default_max_ulps()
+    }
+
+    #[inline]
+    fn relative_eq(&self, other: &Self, epsilon: E, max_relative: E) -> bool {
+        S::Scalar::relative_eq(&self.scale, &other.scale, epsilon, max_relative) &&
+        R::relative_eq(&self.rot, &other.rot, epsilon, max_relative) &&
+        S::relative_eq(&self.disp, &other.disp, epsilon, max_relative)
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: E, max_ulps: u32) -> bool {
+        S::Scalar::ulps_eq(&self.scale, &other.scale, epsilon, max_ulps) &&
+        R::ulps_eq(&self.rot, &other.rot, epsilon, max_ulps) &&
+        S::ulps_eq(&self.disp, &other.disp, epsilon, max_ulps)
     }
 }
 
@@ -174,41 +197,11 @@ mod eders_ser {
         fn serialize<S>(&self, serializer: &mut S) -> Result<(), S::Error>
             where S: serde::Serializer
         {
-            serializer.serialize_struct("Decomposed", DecomposedVisitor {
-                value: self,
-                state: 0,
-            })
-        }
-    }
-
-    struct DecomposedVisitor<'a, V: 'a + VectorSpace, R: 'a> {
-        value: &'a Decomposed<V, R>,
-        state: u8,
-    }
-
-    impl<'a, V: 'a + VectorSpace, R> serde::ser::MapVisitor for DecomposedVisitor<'a, V, R>
-        where V: Serialize, V::Scalar: Serialize, R: Serialize
-    {
-        fn visit<S>(&mut self, serializer: &mut S) -> Result<Option<()>, S::Error>
-            where S: serde::Serializer
-        {
-            match self.state {
-                0 => {
-                    self.state += 1;
-                    Ok(Some(try!(serializer.serialize_struct_elt("scale", &self.value.scale))))
-                },
-                1 => {
-                    self.state += 1;
-                    Ok(Some(try!(serializer.serialize_struct_elt("rot", &self.value.rot))))
-                },
-                2 => {
-                    self.state += 1;
-                    Ok(Some(try!(serializer.serialize_struct_elt("disp", &self.value.disp))))
-                },
-                _ => {
-                    Ok(None)
-                },
-            }
+            let mut state = try!(serializer.serialize_struct("Decomposed", 3));
+            try!(serializer.serialize_struct_elt(&mut state, "scale", &self.scale));
+            try!(serializer.serialize_struct_elt(&mut state, "rot", &self.rot));
+            try!(serializer.serialize_struct_elt(&mut state, "disp", &self.disp));
+            serializer.serialize_struct_end(state)
         }
     }
 }
diff --git a/src/vector.rs b/src/vector.rs
index d3c7997..10369f1 100644
--- a/src/vector.rs
+++ b/src/vector.rs
@@ -172,11 +172,31 @@ macro_rules! impl_vector {
         }
 
         impl<S: BaseFloat> ApproxEq for $VectorN<S> {
-            type Epsilon = S;
+            type Epsilon = S::Epsilon;
 
             #[inline]
-            fn approx_eq_eps(&self, other: &$VectorN<S>, epsilon: &S) -> bool {
-                $(self.$field.approx_eq_eps(&other.$field, epsilon))&&+
+            fn default_epsilon() -> S::Epsilon {
+                S::default_epsilon()
+            }
+
+            #[inline]
+            fn default_max_relative() -> S::Epsilon {
+                S::default_max_relative()
+            }
+
+            #[inline]
+            fn default_max_ulps() -> u32 {
+                S::default_max_ulps()
+            }
+
+            #[inline]
+            fn relative_eq(&self, other: &Self, epsilon: S::Epsilon, max_relative: S::Epsilon) -> bool {
+                $(S::relative_eq(&self.$field, &other.$field, epsilon, max_relative))&&+
+            }
+
+            #[inline]
+            fn ulps_eq(&self, other: &Self, epsilon: S::Epsilon, max_ulps: u32) -> bool {
+                $(S::ulps_eq(&self.$field, &other.$field, epsilon, max_ulps))&&+
             }
         }
 
diff --git a/tests/angle.rs b/tests/angle.rs
index 15b3a0b..02c21bb 100644
--- a/tests/angle.rs
+++ b/tests/angle.rs
@@ -13,26 +13,27 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
+#[macro_use]
+extern crate approx;
 extern crate cgmath;
 
 use cgmath::{Rad, Deg};
-use cgmath::ApproxEq;
 
 #[test]
 fn conv() {
     let angle: Rad<_> = Deg(-5.0f64).into();
     let angle: Deg<_> = angle.into();
-    assert!(angle.approx_eq(&Deg(-5.0f64)));
+    assert_ulps_eq!(&angle, &Deg(-5.0f64));
 
     let angle: Rad<_> = Deg(30.0f64).into();
     let angle: Deg<_> = angle.into();
-    assert!(angle.approx_eq(&Deg(30.0f64)));
+    assert_ulps_eq!(&angle, &Deg(30.0f64));
 
     let angle: Deg<_> = Rad(-5.0f64).into();
     let angle: Rad<_> = angle.into();
-    assert!(angle.approx_eq(&Rad(-5.0f64)));
+    assert_ulps_eq!(&angle, &Rad(-5.0f64));
 
     let angle: Deg<_> = Rad(30.0f64).into();
     let angle: Rad<_> = angle.into();
-    assert!(angle.approx_eq(&Rad(30.0f64)));
+    assert_ulps_eq!(&angle, &Rad(30.0f64));
 }
diff --git a/tests/approx.rs b/tests/approx.rs
deleted file mode 100644
index 2f94a31..0000000
--- a/tests/approx.rs
+++ /dev/null
@@ -1,41 +0,0 @@
-// Copyright 2014 The CGMath Developers. For a full listing of the authors,
-// refer to the Cargo.toml file at the top-level directory of this distribution.
-//
-// Licensed under the Apache License, Version 2.0 (the "License");
-// you may not use this file except in compliance with the License.
-// You may obtain a copy of the License at
-//
-//     http://www.apache.org/licenses/LICENSE-2.0
-//
-// Unless required by applicable law or agreed to in writing, software
-// distributed under the License is distributed on an "AS IS" BASIS,
-// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-// See the License for the specific language governing permissions and
-// limitations under the License.
-
-#[macro_use]
-extern crate cgmath;
-
-use cgmath::*;
-
-#[test]
-fn macro_assert_approx_eq_eps() {
-    assert_approx_eq_eps!(1.0f64, 1.001, 0.01);
-}
-
-#[test]
-#[should_panic]
-fn macro_assert_approx_eq_eps_fail() {
-    assert_approx_eq_eps!(1.0f32, 1.02, 0.01);
-}
-
-#[test]
-fn macro_assert_approx_eq() {
-    assert_approx_eq!(7.0f32 / 5.0, 1.4);
-}
-
-#[test]
-#[should_panic]
-fn macro_assert_approx_eq_fail() {
-    assert_approx_eq!(1.0f64 / 3.0, 0.333);
-}
diff --git a/tests/matrix.rs b/tests/matrix.rs
index abd6b6c..f888f08 100644
--- a/tests/matrix.rs
+++ b/tests/matrix.rs
@@ -13,6 +13,8 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
+#[macro_use]
+extern crate approx;
 #[macro_use]
 extern crate cgmath;
 
@@ -153,15 +155,15 @@ pub mod matrix2 {
     fn test_from_angle() {
         // Rotate the vector (1, 0) by π/2 radians to the vector (0, 1)
         let rot1 = Matrix2::from_angle(Rad(0.5f64 * f64::consts::PI));
-        assert_approx_eq!(rot1 * Vector2::unit_x(), &Vector2::unit_y());
+        assert_ulps_eq!(rot1 * Vector2::unit_x(), &Vector2::unit_y());
 
         // Rotate the vector (-1, 0) by -π/2 radians to the vector (0, 1)
         let rot2 = -rot1;
-        assert_approx_eq!(rot2 * -Vector2::unit_x(), &Vector2::unit_y());
+        assert_ulps_eq!(rot2 * -Vector2::unit_x(), &Vector2::unit_y());
 
         // Rotate the vector (1, 1) by π radians to the vector (-1, -1)
         let rot3: Matrix2<f64> = Matrix2::from_angle(Rad(f64::consts::PI));
-        assert_approx_eq!(rot3 * Vector2::new(1.0, 1.0), &Vector2::new(-1.0, -1.0));
+        assert_ulps_eq!(rot3 * Vector2::new(1.0, 1.0), &Vector2::new(-1.0, -1.0));
     }
 }
 
@@ -319,7 +321,7 @@ pub mod matrix3 {
         fn check_from_axis_angle_x(pitch: Rad<f32>) {
             let found = Matrix3::from_angle_x(pitch);
             let expected = Matrix3::from(Euler { x: pitch, y: Rad(0.0), z: Rad(0.0) });
-            assert_approx_eq_eps!(found, expected, 0.001);
+            assert_relative_eq!(found, expected, epsilon = 0.001);
         }
 
         #[test] fn test_zero()      { check_from_axis_angle_x(Rad(0.0)); }
@@ -333,7 +335,7 @@ pub mod matrix3 {
         fn check_from_axis_angle_y(yaw: Rad<f32>) {
             let found = Matrix3::from_angle_y(yaw);
             let expected = Matrix3::from(Euler { x: Rad(0.0), y: yaw, z: Rad(0.0) });
-            assert_approx_eq_eps!(found, expected, 0.001);
+            assert_relative_eq!(found, expected, epsilon = 0.001);
         }
 
         #[test] fn test_zero()      { check_from_axis_angle_y(Rad(0.0)); }
@@ -347,7 +349,7 @@ pub mod matrix3 {
         fn check_from_axis_angle_z(roll: Rad<f32>) {
             let found = Matrix3::from_angle_z(roll);
             let expected = Matrix3::from(Euler { x: Rad(0.0), y: Rad(0.0), z: roll });
-            assert_approx_eq_eps!(found, expected, 0.001);
+            assert_relative_eq!(found, expected, epsilon = 0.001);
         }
 
         #[test] fn test_zero()      { check_from_axis_angle_z(Rad(0.0)); }
@@ -362,7 +364,7 @@ pub mod matrix3 {
             fn check_from_axis_angle_x(pitch: Rad<f32>) {
                 let found = Matrix3::from_axis_angle(Vector3::unit_x(), pitch);
                 let expected = Matrix3::from(Euler { x: pitch, y: Rad(0.0), z: Rad(0.0) });
-                assert_approx_eq_eps!(found, expected, 0.001);
+                assert_relative_eq!(found, expected, epsilon = 0.001);
             }
 
             #[test] fn test_zero()      { check_from_axis_angle_x(Rad(0.0)); }
@@ -376,7 +378,7 @@ pub mod matrix3 {
             fn check_from_axis_angle_y(yaw: Rad<f32>) {
                 let found = Matrix3::from_axis_angle(Vector3::unit_y(), yaw);
                 let expected = Matrix3::from(Euler { x: Rad(0.0), y: yaw, z: Rad(0.0) });
-                assert_approx_eq_eps!(found, expected, 0.001);
+                assert_relative_eq!(found, expected, epsilon = 0.001);
             }
 
             #[test] fn test_zero()      { check_from_axis_angle_y(Rad(0.0)); }
@@ -390,7 +392,7 @@ pub mod matrix3 {
             fn check_from_axis_angle_z(roll: Rad<f32>) {
                 let found = Matrix3::from_axis_angle(Vector3::unit_z(), roll);
                 let expected = Matrix3::from(Euler { x: Rad(0.0), y: Rad(0.0), z: roll });
-                assert_approx_eq_eps!(found, expected, 0.001);
+                assert_relative_eq!(found, expected, epsilon = 0.001);
             }
 
             #[test] fn test_zero()      { check_from_axis_angle_z(Rad(0.0)); }
@@ -407,10 +409,10 @@ pub mod matrix3 {
             let vec = vec3(0.0, 0.0, 1.0);
 
             let rot = Matrix3::from(Euler::new(Deg(90.0), Deg(0.0), Deg(0.0)));
-            assert_approx_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
 
             let rot = Matrix3::from(Euler::new(Deg(-90.0), Deg(0.0), Deg(0.0)));
-            assert_approx_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
         }
 
         #[test]
@@ -418,10 +420,10 @@ pub mod matrix3 {
             let vec = vec3(0.0, 0.0, 1.0);
 
             let rot = Matrix3::from(Euler::new(Deg(0.0), Deg(90.0), Deg(0.0)));
-            assert_approx_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
 
             let rot = Matrix3::from(Euler::new(Deg(0.0), Deg(-90.0), Deg(0.0)));
-            assert_approx_eq!(vec3(-1.0, 0.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(-1.0, 0.0, 0.0), rot * vec);
         }
 
         #[test]
@@ -429,10 +431,10 @@ pub mod matrix3 {
             let vec = vec3(1.0, 0.0, 0.0);
 
             let rot = Matrix3::from(Euler::new(Deg(0.0), Deg(0.0), Deg(90.0)));
-            assert_approx_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
 
             let rot = Matrix3::from(Euler::new(Deg(0.0), Deg(0.0), Deg(-90.0)));
-            assert_approx_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
         }
 
 
@@ -442,7 +444,7 @@ pub mod matrix3 {
             let vec = vec3(0.0, 1.0, 0.0);
 
             let rot = Matrix3::from(Euler::new(Deg(90.0), Deg(90.0), Deg(0.0)));
-            assert_approx_eq!(vec3(0.0, 0.0, 1.0), rot * vec);
+            assert_ulps_eq!(vec3(0.0, 0.0, 1.0), rot * vec);
         }
 
         // tests that the Z rotation is done after the Y
@@ -451,7 +453,7 @@ pub mod matrix3 {
             let vec = vec3(0.0, 0.0, 1.0);
 
             let rot = Matrix3::from(Euler::new(Deg(0.0), Deg(90.0), Deg(90.0)));
-            assert_approx_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
         }
     }
 
@@ -464,7 +466,7 @@ pub mod matrix3 {
 
             let rot = Matrix3::from_angle_x(Deg(90.0));
             println!("x mat: {:?}", rot);
-            assert_approx_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
         }
 
         #[test]
@@ -472,7 +474,7 @@ pub mod matrix3 {
             let vec = vec3(0.0, 0.0, 1.0);
 
             let rot = Matrix3::from_angle_y(Deg(90.0));
-            assert_approx_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
         }
 
         #[test]
@@ -480,7 +482,7 @@ pub mod matrix3 {
             let vec = vec3(1.0, 0.0, 0.0);
 
             let rot = Matrix3::from_angle_z(Deg(90.0));
-            assert_approx_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
         }
 
         #[test]
@@ -488,7 +490,7 @@ pub mod matrix3 {
             let vec = vec3(0.0, 0.0, 1.0);
 
             let rot = Matrix3::from_axis_angle(vec3(1.0, 1.0, 0.0).normalize(), Deg(90.0));
-            assert_approx_eq!(vec3(2.0f32.sqrt() / 2.0, -2.0f32.sqrt() / 2.0, 0.0), rot * vec);
+            assert_ulps_eq!(vec3(2.0f32.sqrt() / 2.0, -2.0f32.sqrt() / 2.0, 0.0), rot * vec);
         }
 
         #[test]
@@ -496,7 +498,7 @@ pub mod matrix3 {
             let vec = vec3(1.0, 0.0, 0.0);
 
             let rot = Matrix3::from_axis_angle(vec3(0.0, 1.0, 1.0).normalize(), Deg(-90.0));
-            assert_approx_eq!(vec3(0.0, -2.0f32.sqrt() / 2.0, 2.0f32.sqrt() / 2.0), rot * vec);
+            assert_ulps_eq!(vec3(0.0, -2.0f32.sqrt() / 2.0, 2.0f32.sqrt() / 2.0), rot * vec);
         }
 
         #[test]
@@ -504,7 +506,7 @@ pub mod matrix3 {
             let vec = vec3(0.0, 1.0, 0.0);
 
             let rot = Matrix3::from_axis_angle(vec3(1.0, 0.0, 1.0).normalize(), Deg(90.0));
-            assert_approx_eq!(vec3(-2.0f32.sqrt() / 2.0, 0.0, 2.0f32.sqrt() / 2.0), rot * vec);
+            assert_ulps_eq!(vec3(-2.0f32.sqrt() / 2.0, 0.0, 2.0f32.sqrt() / 2.0), rot * vec);
         }
     }
 }
@@ -643,11 +645,11 @@ pub mod matrix4 {
     fn test_invert() {
         assert!(Matrix4::<f64>::identity().invert().unwrap().is_identity());
 
-        assert!(C.invert().unwrap().approx_eq(&(
+        assert_ulps_eq!(&C.invert().unwrap(), &(
                 Matrix4::new( 5.0f64, -4.0f64,  1.0f64,  0.0f64,
                              -4.0f64,  8.0f64, -4.0f64,  0.0f64,
                               4.0f64, -8.0f64,  4.0f64,  8.0f64,
-                             -3.0f64,  4.0f64,  1.0f64, -8.0f64) * 0.125f64)));
+                             -3.0f64,  4.0f64,  1.0f64, -8.0f64) * 0.125f64));
 
         let mat_c = Matrix4::new(-0.131917f64, -0.76871f64,   0.625846f64, 0.0f64,
                                  -0.,        0.631364f64,  0.775487f64, 0.0f64,
@@ -714,7 +716,7 @@ pub mod matrix4 {
             let matrix_long = Matrix3::from(quaternion);
             let matrix_long = Matrix4::from(matrix_long);
 
-            assert_approx_eq!(matrix_short, matrix_long);
+            assert_ulps_eq!(matrix_short, matrix_long);
         }
     }
 }
diff --git a/tests/point.rs b/tests/point.rs
index fcd2714..dc896d7 100644
--- a/tests/point.rs
+++ b/tests/point.rs
@@ -13,11 +13,11 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
-
+#[macro_use]
+extern crate approx;
 extern crate cgmath;
 
 use cgmath::{Point2, Point3};
-use cgmath::ApproxEq;
 
 macro_rules! impl_test_mul {
     ($PointN:ident { $($field:ident),+ }, $s:expr, $v:expr) => (
@@ -54,7 +54,7 @@ macro_rules! impl_test_rem {
 #[test]
 fn test_homogeneous() {
 	let p = Point3::new(1.0f64, 2.0f64, 3.0f64);
-    assert!(p.approx_eq(&Point3::from_homogeneous(p.to_homogeneous())));
+    assert_ulps_eq!(&p, &Point3::from_homogeneous(p.to_homogeneous()));
 }
 
 #[test]
@@ -74,4 +74,3 @@ fn test_rem() {
     impl_test_rem!(Point3 { x, y, z }, 2.0f32, Point3::new(2.0f32, 4.0, 6.0));
     impl_test_rem!(Point2 { x, y }, 2.0f32, Point2::new(2.0f32, 4.0));
 }
-
diff --git a/tests/quaternion.rs b/tests/quaternion.rs
index a99ddda..f59911d 100644
--- a/tests/quaternion.rs
+++ b/tests/quaternion.rs
@@ -13,6 +13,8 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
+#[macro_use]
+extern crate approx;
 #[macro_use]
 extern crate cgmath;
 
@@ -58,7 +60,7 @@ mod to_from_euler {
     use cgmath::*;
 
     fn check_euler(rotation: Euler<Rad<f32>>) {
-        assert_approx_eq_eps!(Euler::from(Quaternion::from(rotation)), rotation, 0.001);
+        assert_relative_eq!(Euler::from(Quaternion::from(rotation)), rotation, epsilon = 0.001);
     }
 
     const HPI: f32 = f32::consts::FRAC_PI_2;
@@ -84,7 +86,7 @@ mod from {
             let matrix3 = Matrix3::from(Euler { x: x, y: y, z: z });
             let quaternion = Quaternion::from(matrix3);
             let quaternion_matrix3 = Matrix3::from(quaternion);
-            assert_approx_eq!(matrix3, quaternion_matrix3);
+            assert_ulps_eq!(matrix3, quaternion_matrix3);
         }
 
         // triggers: trace >= S::zero()
@@ -120,7 +122,7 @@ mod arc {
     fn test(src: Vector3<f32>, dst: Vector3<f32>) {
         let q = Quaternion::from_arc(src, dst, None);
         let v = q.rotate_vector(src);
-        assert_approx_eq!(v.normalize(), dst.normalize());
+        assert_ulps_eq!(v.normalize(), dst.normalize());
     }
 
     #[test]
@@ -145,7 +147,7 @@ mod arc {
     fn test_ortho() {
         let q: Quaternion<f32> = Quaternion::from_arc(Vector3::unit_x(), Vector3::unit_y(), None);
         let q2 = Quaternion::from_axis_angle(Vector3::unit_z(), Rad::turn_div_4());
-        assert_approx_eq!(q, q2);
+        assert_ulps_eq!(q, q2);
     }
 }
 
@@ -157,10 +159,10 @@ mod rotate_from_euler {
         let vec = vec3(0.0, 0.0, 1.0);
 
         let rot = Quaternion::from(Euler::new(Deg(90.0), Deg(0.0), Deg(0.0)));
-        assert_approx_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
 
         let rot = Quaternion::from(Euler::new(Deg(-90.0), Deg(0.0), Deg(0.0)));
-        assert_approx_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
     }
 
     #[test]
@@ -168,10 +170,10 @@ mod rotate_from_euler {
         let vec = vec3(0.0, 0.0, 1.0);
 
         let rot = Quaternion::from(Euler::new(Deg(0.0), Deg(90.0), Deg(0.0)));
-        assert_approx_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
 
         let rot = Quaternion::from(Euler::new(Deg(0.0), Deg(-90.0), Deg(0.0)));
-        assert_approx_eq!(vec3(-1.0, 0.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(-1.0, 0.0, 0.0), rot * vec);
     }
 
     #[test]
@@ -179,10 +181,10 @@ mod rotate_from_euler {
         let vec = vec3(1.0, 0.0, 0.0);
 
         let rot = Quaternion::from(Euler::new(Deg(0.0), Deg(0.0), Deg(90.0)));
-        assert_approx_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
 
         let rot = Quaternion::from(Euler::new(Deg(0.0), Deg(0.0), Deg(-90.0)));
-        assert_approx_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
     }
 
 
@@ -192,7 +194,7 @@ mod rotate_from_euler {
         let vec = vec3(0.0, 1.0, 0.0);
 
         let rot = Quaternion::from(Euler::new(Deg(90.0), Deg(90.0), Deg(0.0)));
-        assert_approx_eq!(vec3(0.0, 0.0, 1.0), rot * vec);
+        assert_ulps_eq!(vec3(0.0, 0.0, 1.0), rot * vec);
     }
 
     // tests that the Z rotation is done after the Y
@@ -201,7 +203,7 @@ mod rotate_from_euler {
         let vec = vec3(0.0, 0.0, 1.0);
 
         let rot = Quaternion::from(Euler::new(Deg(0.0), Deg(90.0), Deg(90.0)));
-        assert_approx_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
     }
 }
 
@@ -213,7 +215,7 @@ mod rotate_from_axis_angle {
         let vec = vec3(0.0, 0.0, 1.0);
 
         let rot = Quaternion::from_angle_x(Deg(90.0));
-        assert_approx_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(0.0, -1.0, 0.0), rot * vec);
     }
 
     #[test]
@@ -221,7 +223,7 @@ mod rotate_from_axis_angle {
         let vec = vec3(0.0, 0.0, 1.0);
 
         let rot = Quaternion::from_angle_y(Deg(90.0));
-        assert_approx_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(1.0, 0.0, 0.0), rot * vec);
     }
 
     #[test]
@@ -229,7 +231,7 @@ mod rotate_from_axis_angle {
         let vec = vec3(1.0, 0.0, 0.0);
 
         let rot = Quaternion::from_angle_z(Deg(90.0));
-        assert_approx_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(0.0, 1.0, 0.0), rot * vec);
     }
 
     #[test]
@@ -237,7 +239,7 @@ mod rotate_from_axis_angle {
         let vec = vec3(0.0, 0.0, 1.0);
 
         let rot = Quaternion::from_axis_angle(vec3(1.0, 1.0, 0.0).normalize(), Deg(90.0));
-        assert_approx_eq!(vec3(2.0f32.sqrt() / 2.0, -2.0f32.sqrt() / 2.0, 0.0), rot * vec);
+        assert_ulps_eq!(vec3(2.0f32.sqrt() / 2.0, -2.0f32.sqrt() / 2.0, 0.0), rot * vec);
     }
 
     #[test]
@@ -245,7 +247,7 @@ mod rotate_from_axis_angle {
         let vec = vec3(1.0, 0.0, 0.0);
 
         let rot = Quaternion::from_axis_angle(vec3(0.0, 1.0, 1.0).normalize(), Deg(-90.0));
-        assert_approx_eq!(vec3(0.0, -2.0f32.sqrt() / 2.0, 2.0f32.sqrt() / 2.0), rot * vec);
+        assert_ulps_eq!(vec3(0.0, -2.0f32.sqrt() / 2.0, 2.0f32.sqrt() / 2.0), rot * vec);
     }
 
     #[test]
@@ -253,6 +255,6 @@ mod rotate_from_axis_angle {
         let vec = vec3(0.0, 1.0, 0.0);
 
         let rot = Quaternion::from_axis_angle(vec3(1.0, 0.0, 1.0).normalize(), Deg(90.0));
-        assert_approx_eq!(vec3(-2.0f32.sqrt() / 2.0, 0.0, 2.0f32.sqrt() / 2.0), rot * vec);
+        assert_ulps_eq!(vec3(-2.0f32.sqrt() / 2.0, 0.0, 2.0f32.sqrt() / 2.0), rot * vec);
     }
 }
diff --git a/tests/transform.rs b/tests/transform.rs
index 1688c69..0c1cfbf 100644
--- a/tests/transform.rs
+++ b/tests/transform.rs
@@ -13,7 +13,8 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
-
+#[macro_use]
+extern crate approx;
 extern crate cgmath;
 
 #[cfg(feature = "eders")]
@@ -31,7 +32,7 @@ fn test_invert() {
     };
     let ti = t.inverse_transform().expect("Expected successful inversion");
     let vt = t.transform_vector(v);
-    assert!(v.approx_eq(&ti.transform_vector(vt)));
+    assert_ulps_eq!(&v, &ti.transform_vector(vt));
 }
 
 #[test]
@@ -42,7 +43,7 @@ fn test_look_at() {
     let t: Decomposed<Vector3<f64>, Quaternion<f64>> = Transform::look_at(eye, center, up);
     let point = Point3::new(1.0f64, 0.0, 0.0);
     let view_point = Point3::new(0.0f64, 1.0, 5.0);
-    assert!(t.transform_point(point).approx_eq(&view_point));
+    assert_ulps_eq!(&t.transform_point(point), &view_point);
 }
 
 #[cfg(feature = "eders")]
@@ -57,5 +58,5 @@ fn test_serialize() {
     let serialized = serde_json::to_string(&t).unwrap();
     let deserialized: Decomposed<Vector3<f64>, Quaternion<f64>> = serde_json::from_str(&serialized).unwrap();
 
-    assert!(t.approx_eq(&deserialized));
+    assert_ulps_eq!(&t, &deserialized);
 }
diff --git a/tests/vector.rs b/tests/vector.rs
index 57a9f9d..307a5de 100644
--- a/tests/vector.rs
+++ b/tests/vector.rs
@@ -13,6 +13,8 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
+#[macro_use]
+extern crate approx;
 #[macro_use]
 extern crate cgmath;
 
@@ -229,30 +231,30 @@ mod test_magnitude {
 
 #[test]
 fn test_angle() {
-    assert!(Vector2::new(1.0f64, 0.0f64).angle(Vector2::new(0.0f64, 1.0f64)).approx_eq( &Rad(f64::consts::FRAC_PI_2) ));
-    assert!(Vector2::new(10.0f64, 0.0f64).angle(Vector2::new(0.0f64, 5.0f64)).approx_eq( &Rad(f64::consts::FRAC_PI_2) ));
-    assert!(Vector2::new(-1.0f64, 0.0f64).angle(Vector2::new(0.0f64, 1.0f64)).approx_eq( &-Rad(f64::consts::FRAC_PI_2) ));
+    assert_ulps_eq!(Vector2::new(1.0f64, 0.0f64).angle(Vector2::new(0.0f64, 1.0f64)), &Rad(f64::consts::FRAC_PI_2));
+    assert_ulps_eq!(Vector2::new(10.0f64, 0.0f64).angle(Vector2::new(0.0f64, 5.0f64)), &Rad(f64::consts::FRAC_PI_2));
+    assert_ulps_eq!(Vector2::new(-1.0f64, 0.0f64).angle(Vector2::new(0.0f64, 1.0f64)), &-Rad(f64::consts::FRAC_PI_2));
 
-    assert!(Vector3::new(1.0f64, 0.0f64, 1.0f64).angle(Vector3::new(1.0f64, 1.0f64, 0.0f64)).approx_eq( &Rad(f64::consts::FRAC_PI_3) ));
-    assert!(Vector3::new(10.0f64, 0.0f64, 10.0f64).angle(Vector3::new(5.0f64, 5.0f64, 0.0f64)).approx_eq( &Rad(f64::consts::FRAC_PI_3) ));
-    assert!(Vector3::new(-1.0f64, 0.0f64, -1.0f64).angle(Vector3::new(1.0f64, -1.0f64, 0.0f64)).approx_eq( &Rad(2.0f64 * f64::consts::FRAC_PI_3) ));
+    assert_ulps_eq!(Vector3::new(1.0f64, 0.0f64, 1.0f64).angle(Vector3::new(1.0f64, 1.0f64, 0.0f64)), &Rad(f64::consts::FRAC_PI_3));
+    assert_ulps_eq!(Vector3::new(10.0f64, 0.0f64, 10.0f64).angle(Vector3::new(5.0f64, 5.0f64, 0.0f64)), &Rad(f64::consts::FRAC_PI_3));
+    assert_ulps_eq!(Vector3::new(-1.0f64, 0.0f64, -1.0f64).angle(Vector3::new(1.0f64, -1.0f64, 0.0f64)), &Rad(2.0f64 * f64::consts::FRAC_PI_3));
 
-    assert!(Vector4::new(1.0f64, 0.0f64, 1.0f64, 0.0f64).angle(Vector4::new(0.0f64, 1.0f64, 0.0f64, 1.0f64)).approx_eq( &Rad(f64::consts::FRAC_PI_2) ));
-    assert!(Vector4::new(10.0f64, 0.0f64, 10.0f64, 0.0f64).angle(Vector4::new(0.0f64, 5.0f64, 0.0f64, 5.0f64)).approx_eq( &Rad(f64::consts::FRAC_PI_2) ));
-    assert!(Vector4::new(-1.0f64, 0.0f64, -1.0f64, 0.0f64).angle(Vector4::new(0.0f64, 1.0f64, 0.0f64, 1.0f64)).approx_eq( &Rad(f64::consts::FRAC_PI_2) ));
+    assert_ulps_eq!(Vector4::new(1.0f64, 0.0f64, 1.0f64, 0.0f64).angle(Vector4::new(0.0f64, 1.0f64, 0.0f64, 1.0f64)), &Rad(f64::consts::FRAC_PI_2));
+    assert_ulps_eq!(Vector4::new(10.0f64, 0.0f64, 10.0f64, 0.0f64).angle(Vector4::new(0.0f64, 5.0f64, 0.0f64, 5.0f64)), &Rad(f64::consts::FRAC_PI_2));
+    assert_ulps_eq!(Vector4::new(-1.0f64, 0.0f64, -1.0f64, 0.0f64).angle(Vector4::new(0.0f64, 1.0f64, 0.0f64, 1.0f64)), &Rad(f64::consts::FRAC_PI_2));
 }
 
 #[test]
 fn test_normalize() {
     // TODO: test normalize_to, normalize_sel.0, and normalize_self_to
-    assert!(Vector2::new(3.0f64, 4.0f64).normalize().approx_eq( &Vector2::new(3.0/5.0, 4.0/5.0) ));
-    assert!(Vector3::new(2.0f64, 3.0f64, 6.0f64).normalize().approx_eq( &Vector3::new(2.0/7.0, 3.0/7.0, 6.0/7.0) ));
-    assert!(Vector4::new(1.0f64, 2.0f64, 4.0f64, 10.0f64).normalize().approx_eq( &Vector4::new(1.0/11.0, 2.0/11.0, 4.0/11.0, 10.0/11.0) ));
+    assert_ulps_eq!(Vector2::new(3.0f64, 4.0f64).normalize(), &Vector2::new(3.0/5.0, 4.0/5.0));
+    assert_ulps_eq!(Vector3::new(2.0f64, 3.0f64, 6.0f64).normalize(), &Vector3::new(2.0/7.0, 3.0/7.0, 6.0/7.0));
+    assert_ulps_eq!(Vector4::new(1.0f64, 2.0f64, 4.0f64, 10.0f64).normalize(), &Vector4::new(1.0/11.0, 2.0/11.0, 4.0/11.0, 10.0/11.0));
 }
 
 #[test]
 fn test_cast() {
-    assert_approx_eq!(Vector2::new(0.9f64, 1.5).cast(), Vector2::new(0.9f32, 1.5));
-    assert_approx_eq!(Vector3::new(1.0f64, 2.4, -3.13).cast(), Vector3::new(1.0f32, 2.4, -3.13));
-    assert_approx_eq!(Vector4::new(13.5f64, -4.6, -8.3, 2.41).cast(), Vector4::new(13.5f32, -4.6, -8.3, 2.41));
+    assert_ulps_eq!(Vector2::new(0.9f64, 1.5).cast(), Vector2::new(0.9f32, 1.5));
+    assert_ulps_eq!(Vector3::new(1.0f64, 2.4, -3.13).cast(), Vector3::new(1.0f32, 2.4, -3.13));
+    assert_ulps_eq!(Vector4::new(13.5f64, -4.6, -8.3, 2.41).cast(), Vector4::new(13.5f32, -4.6, -8.3, 2.41));
 }