diff --git a/tests/core/blas.rs b/tests/core/blas.rs
index 57fab2c6..82239af0 100644
--- a/tests/core/blas.rs
+++ b/tests/core/blas.rs
@@ -1,9 +1,8 @@
-use std::cmp;
+#![cfg(feature = "arbitrary")]
 
+use std::cmp;
 use na::{DVector, DMatrix};
 
-
-#[cfg(feature = "arbitrary")]
 quickcheck! {
     /*
      *
diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs
index 2ebab858..996d1d64 100644
--- a/tests/core/matrix.rs
+++ b/tests/core/matrix.rs
@@ -1,8 +1,5 @@
 use num::{Zero, One};
 use num::Float;
-use std::fmt::Display;
-
-use alga::linear::FiniteDimInnerSpace;
 
 use na::{self,
          DVector, DMatrix,
@@ -806,93 +803,100 @@ mod normalization_tests {
     }
 }
 
+#[cfg(feature = "arbitrary")]
 // FIXME: move this to alga ?
-macro_rules! finite_dim_inner_space_test(
-    ($($Vector: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$(
-        #[cfg(feature = "arbitrary")]
-        quickcheck!{
-            fn $orthonormal_subspace(vs: Vec<$Vector<f64>>) -> bool {
-                let mut given_basis = vs.clone();
-                let given_basis_dim = $Vector::orthonormalize(&mut given_basis[..]);
-                let mut ortho_basis = Vec::new();
-                $Vector::orthonormal_subspace_basis(
-                    &given_basis[.. given_basis_dim],
-                    |e| { ortho_basis.push(*e); true }
-                );
+mod finite_dim_inner_space_tests {
+    use super::*;
+    use std::fmt::Display;
+    use alga::linear::FiniteDimInnerSpace;
 
-                if !is_subspace_basis(&ortho_basis[..]) {
-                    return false;
+    macro_rules! finite_dim_inner_space_test(
+        ($($Vector: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$(
+            quickcheck!{
+                fn $orthonormal_subspace(vs: Vec<$Vector<f64>>) -> bool {
+                    let mut given_basis = vs.clone();
+                    let given_basis_dim = $Vector::orthonormalize(&mut given_basis[..]);
+                    let mut ortho_basis = Vec::new();
+                    $Vector::orthonormal_subspace_basis(
+                        &given_basis[.. given_basis_dim],
+                        |e| { ortho_basis.push(*e); true }
+                    );
+
+                    if !is_subspace_basis(&ortho_basis[..]) {
+                        return false;
+                    }
+
+                    for v in vs {
+                        for b in &ortho_basis {
+                            if !relative_eq!(v.dot(b), 0.0, epsilon = 1.0e-7) {
+                                println!("Found dot product: {} · {} = {}", v, b, v.dot(b));
+                                return false;
+                            }
+                        }
+                    }
+
+                    true
                 }
 
-                for v in vs {
-                    for b in &ortho_basis {
-                        if !relative_eq!(v.dot(b), 0.0, epsilon = 1.0e-7) {
-                            println!("Found dot product: {} · {} = {}", v, b, v.dot(b));
+                fn $orthonormalization(vs: Vec<$Vector<f64>>) -> bool {
+                    let mut basis = vs.clone();
+                    let subdim = $Vector::orthonormalize(&mut basis[..]);
+
+                    if !is_subspace_basis(&basis[.. subdim]) {
+                        return false;
+                    }
+
+                    for mut e in vs {
+                        for b in &basis[.. subdim] {
+                            e -= e.dot(b) * b
+                        }
+
+                        // Any element of `e` must be a linear combination of the basis elements.
+                        if !relative_eq!(e.norm(), 0.0, epsilon = 1.0e-7) {
+                            println!("Orthonormalization; element decomposition failure: {}", e);
+                            println!("... the non-zero norm is: {}", e.norm());
                             return false;
                         }
                     }
-                }
 
-                true
+                    true
+                }
+            }
+        )*}
+    );
+
+    finite_dim_inner_space_test!(
+        Vector1, orthonormal_subspace_basis1, orthonormalize1;
+        Vector2, orthonormal_subspace_basis2, orthonormalize2;
+        Vector3, orthonormal_subspace_basis3, orthonormalize3;
+        Vector4, orthonormal_subspace_basis4, orthonormalize4;
+        Vector5, orthonormal_subspace_basis5, orthonormalize5;
+        Vector6, orthonormal_subspace_basis6, orthonormalize6;
+    );
+
+    /*
+     *
+     * Helper functions.
+     *
+     */
+    #[cfg(feature = "arbitrary")]
+    fn is_subspace_basis<T: FiniteDimInnerSpace<Real = f64> + Display>(vs: &[T]) -> bool {
+        for i in 0 .. vs.len() {
+            // Basis elements must be normalized.
+            if !relative_eq!(vs[i].norm(), 1.0, epsilon = 1.0e-7) {
+                println!("Non-zero basis element norm: {}", vs[i].norm());
+                return false;
             }
 
-            fn $orthonormalization(vs: Vec<$Vector<f64>>) -> bool {
-                let mut basis = vs.clone();
-                let subdim = $Vector::orthonormalize(&mut basis[..]);
-
-                if !is_subspace_basis(&basis[.. subdim]) {
-                    return false;
+            for j in 0 .. i {
+                // Basis elements must be orthogonal.
+                if !relative_eq!(vs[i].dot(&vs[j]), 0.0, epsilon = 1.0e-7) {
+                    println!("Non-orthogonal basis elements: {} · {} = {}", vs[i], vs[j], vs[i].dot(&vs[j]));
+                    return false
                 }
-
-                for mut e in vs {
-                    for b in &basis[.. subdim] {
-                        e -= e.dot(b) * b
-                    }
-
-                    // Any element of `e` must be a linear combination of the basis elements.
-                    if !relative_eq!(e.norm(), 0.0, epsilon = 1.0e-7) {
-                        println!("Orthonormalization; element decomposition failure: {}", e);
-                        println!("... the non-zero norm is: {}", e.norm());
-                        return false;
-                    }
-                }
-
-                true
             }
         }
-    )*}
-);
 
-finite_dim_inner_space_test!(
-    Vector1, orthonormal_subspace_basis1, orthonormalize1;
-    Vector2, orthonormal_subspace_basis2, orthonormalize2;
-    Vector3, orthonormal_subspace_basis3, orthonormalize3;
-    Vector4, orthonormal_subspace_basis4, orthonormalize4;
-    Vector5, orthonormal_subspace_basis5, orthonormalize5;
-    Vector6, orthonormal_subspace_basis6, orthonormalize6;
-);
-
-/*
- *
- * Helper functions.
- *
- */
-fn is_subspace_basis<T: FiniteDimInnerSpace<Real = f64> + Display>(vs: &[T]) -> bool {
-    for i in 0 .. vs.len() {
-        // Basis elements must be normalized.
-        if !relative_eq!(vs[i].norm(), 1.0, epsilon = 1.0e-7) {
-            println!("Non-zero basis element norm: {}", vs[i].norm());
-            return false;
-        }
-
-        for j in 0 .. i {
-            // Basis elements must be orthogonal.
-            if !relative_eq!(vs[i].dot(&vs[j]), 0.0, epsilon = 1.0e-7) {
-                println!("Non-orthogonal basis elements: {} · {} = {}", vs[i], vs[j], vs[i].dot(&vs[j]));
-                return false
-            }
-        }
+        true
     }
-
-    true
 }
diff --git a/tests/geometry/projection.rs b/tests/geometry/projection.rs
index 84f3cf70..7be0d49d 100644
--- a/tests/geometry/projection.rs
+++ b/tests/geometry/projection.rs
@@ -1,4 +1,4 @@
-use na::{Point3, Perspective3, Orthographic3};
+use na::{Perspective3, Orthographic3};
 
 #[test]
 fn perspective_inverse() {
@@ -22,22 +22,26 @@ fn orthographic_inverse() {
 
 
 #[cfg(feature = "arbitrary")]
-quickcheck!{
-    fn perspective_project_unproject(pt: Point3<f64>) -> bool {
-        let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0);
+mod quickcheck_tests {
+    use na::{Point3, Perspective3, Orthographic3};
 
-        let projected   = proj.project_point(&pt);
-        let unprojected = proj.unproject_point(&projected);
+    quickcheck!{
+        fn perspective_project_unproject(pt: Point3<f64>) -> bool {
+            let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0);
 
-        relative_eq!(pt, unprojected, epsilon = 1.0e-7)
-    }
+            let projected   = proj.project_point(&pt);
+            let unprojected = proj.unproject_point(&projected);
 
-    fn orthographic_project_unproject(pt: Point3<f64>) -> bool {
-        let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0);
+            relative_eq!(pt, unprojected, epsilon = 1.0e-7)
+        }
 
-        let projected   = proj.project_point(&pt);
-        let unprojected = proj.unproject_point(&projected);
+        fn orthographic_project_unproject(pt: Point3<f64>) -> bool {
+            let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0);
 
-        relative_eq!(pt, unprojected, epsilon = 1.0e-7)
+            let projected   = proj.project_point(&pt);
+            let unprojected = proj.unproject_point(&projected);
+
+            relative_eq!(pt, unprojected, epsilon = 1.0e-7)
+        }
     }
 }
diff --git a/tests/geometry/rotation.rs b/tests/geometry/rotation.rs
index a6e6daf6..a2e91dec 100644
--- a/tests/geometry/rotation.rs
+++ b/tests/geometry/rotation.rs
@@ -1,6 +1,4 @@
-use std::f64;
-use alga::general::Real;
-use na::{self, Vector2, Vector3, Rotation2, Rotation3, Unit};
+use na::{Vector2, Vector3};
 
 #[test]
 fn angle_2() {
@@ -19,187 +17,193 @@ fn angle_3() {
 }
 
 #[cfg(feature = "arbitrary")]
-quickcheck!(
-    /*
-     *
-     * Euler angles.
-     *
-     */
-    fn from_euler_angles(r: f64, p: f64, y: f64) -> bool {
-        let roll  = Rotation3::from_euler_angles(r, 0.0, 0.0);
-        let pitch = Rotation3::from_euler_angles(0.0, p, 0.0);
-        let yaw   = Rotation3::from_euler_angles(0.0, 0.0, y);
+mod quickcheck_tests {
+    use std::f64;
+    use alga::general::Real;
+    use na::{self, Vector2, Vector3, Rotation2, Rotation3, Unit};
 
-        let rpy = Rotation3::from_euler_angles(r, p, y);
+    quickcheck! {
+        /*
+         *
+         * Euler angles.
+         *
+         */
+        fn from_euler_angles(r: f64, p: f64, y: f64) -> bool {
+            let roll  = Rotation3::from_euler_angles(r, 0.0, 0.0);
+            let pitch = Rotation3::from_euler_angles(0.0, p, 0.0);
+            let yaw   = Rotation3::from_euler_angles(0.0, 0.0, y);
 
-        roll[(0, 0)]  == 1.0 && // rotation wrt. x axis.
-        pitch[(1, 1)] == 1.0 && // rotation wrt. y axis.
-        yaw[(2, 2)]   == 1.0 && // rotation wrt. z axis.
-        yaw * pitch * roll == rpy
-    }
+            let rpy = Rotation3::from_euler_angles(r, p, y);
 
-    fn to_euler_angles(r: f64, p: f64, y: f64) -> bool {
-        let rpy = Rotation3::from_euler_angles(r, p, y);
-        let (roll, pitch, yaw) = rpy.to_euler_angles();
-        relative_eq!(Rotation3::from_euler_angles(roll, pitch, yaw), rpy, epsilon = 1.0e-7)
-    }
-
-    fn to_euler_angles_gimble_lock(r: f64, y: f64) -> bool {
-        let pos = Rotation3::from_euler_angles(r,  f64::frac_pi_2(), y);
-        let neg = Rotation3::from_euler_angles(r, -f64::frac_pi_2(), y);
-        let (pos_r, pos_p, pos_y) = pos.to_euler_angles();
-        let (neg_r, neg_p, neg_y) = neg.to_euler_angles();
-        relative_eq!(Rotation3::from_euler_angles(pos_r, pos_p, pos_y), pos, epsilon = 1.0e-7) &&
-        relative_eq!(Rotation3::from_euler_angles(neg_r, neg_p, neg_y), neg, epsilon = 1.0e-7)
-    }
-
-    /*
-     *
-     * Inversion is transposition.
-     *
-     */
-    fn rotation_inv_3(a: Rotation3<f64>) -> bool {
-        let ta = a.transpose();
-        let ia = a.inverse();
-
-        ta == ia &&
-        relative_eq!(&ta * &a,  Rotation3::identity(), epsilon = 1.0e-7) &&
-        relative_eq!(&ia *  a,  Rotation3::identity(), epsilon = 1.0e-7) &&
-        relative_eq!(  a * &ta, Rotation3::identity(), epsilon = 1.0e-7) &&
-        relative_eq!(  a *  ia, Rotation3::identity(), epsilon = 1.0e-7)
-    }
-
-    fn rotation_inv_2(a: Rotation2<f64>) -> bool {
-        let ta = a.transpose();
-        let ia = a.inverse();
-
-        ta == ia &&
-        relative_eq!(&ta * &a,  Rotation2::identity(), epsilon = 1.0e-7) &&
-        relative_eq!(&ia *  a,  Rotation2::identity(), epsilon = 1.0e-7) &&
-        relative_eq!(  a * &ta, Rotation2::identity(), epsilon = 1.0e-7) &&
-        relative_eq!(  a *  ia, Rotation2::identity(), epsilon = 1.0e-7)
-    }
-
-    /*
-     *
-     * Angle between vectors.
-     *
-     */
-    fn angle_is_commutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
-        a.angle(&b) == b.angle(&a)
-    }
-
-    fn angle_is_commutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
-        a.angle(&b) == b.angle(&a)
-    }
-
-    /*
-     *
-     * Rotation matrix between vectors.
-     *
-     */
-    fn rotation_between_is_anticommutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
-        let rab = Rotation2::rotation_between(&a, &b);
-        let rba = Rotation2::rotation_between(&b, &a);
-
-        relative_eq!(rab * rba, Rotation2::identity())
-    }
-
-    fn rotation_between_is_anticommutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
-        let rots = (Rotation3::rotation_between(&a, &b), Rotation3::rotation_between(&b, &a));
-        if let (Some(rab), Some(rba)) = rots {
-            relative_eq!(rab * rba, Rotation3::identity(), epsilon = 1.0e-7)
+            roll[(0, 0)]  == 1.0 && // rotation wrt. x axis.
+            pitch[(1, 1)] == 1.0 && // rotation wrt. y axis.
+            yaw[(2, 2)]   == 1.0 && // rotation wrt. z axis.
+            yaw * pitch * roll == rpy
         }
-        else {
-            true
+
+        fn to_euler_angles(r: f64, p: f64, y: f64) -> bool {
+            let rpy = Rotation3::from_euler_angles(r, p, y);
+            let (roll, pitch, yaw) = rpy.to_euler_angles();
+            relative_eq!(Rotation3::from_euler_angles(roll, pitch, yaw), rpy, epsilon = 1.0e-7)
         }
-    }
 
-    fn rotation_between_is_identity(v2: Vector2<f64>, v3: Vector3<f64>) -> bool {
-        let vv2 = 3.42 * v2;
-        let vv3 = 4.23 * v3;
-
-        relative_eq!(v2.angle(&vv2), 0.0, epsilon = 1.0e-7) &&
-        relative_eq!(v3.angle(&vv3), 0.0, epsilon = 1.0e-7) &&
-        relative_eq!(Rotation2::rotation_between(&v2, &vv2), Rotation2::identity()) &&
-        Rotation3::rotation_between(&v3, &vv3).unwrap() == Rotation3::identity()
-    }
-
-    fn rotation_between_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
-        if !relative_eq!(a.angle(&b), 0.0, epsilon = 1.0e-7) {
-            let r = Rotation2::rotation_between(&a, &b);
-            relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7)
+        fn to_euler_angles_gimble_lock(r: f64, y: f64) -> bool {
+            let pos = Rotation3::from_euler_angles(r,  f64::frac_pi_2(), y);
+            let neg = Rotation3::from_euler_angles(r, -f64::frac_pi_2(), y);
+            let (pos_r, pos_p, pos_y) = pos.to_euler_angles();
+            let (neg_r, neg_p, neg_y) = neg.to_euler_angles();
+            relative_eq!(Rotation3::from_euler_angles(pos_r, pos_p, pos_y), pos, epsilon = 1.0e-7) &&
+            relative_eq!(Rotation3::from_euler_angles(neg_r, neg_p, neg_y), neg, epsilon = 1.0e-7)
         }
-        else {
-            true
+
+        /*
+         *
+         * Inversion is transposition.
+         *
+         */
+        fn rotation_inv_3(a: Rotation3<f64>) -> bool {
+            let ta = a.transpose();
+            let ia = a.inverse();
+
+            ta == ia &&
+            relative_eq!(&ta * &a,  Rotation3::identity(), epsilon = 1.0e-7) &&
+            relative_eq!(&ia *  a,  Rotation3::identity(), epsilon = 1.0e-7) &&
+            relative_eq!(  a * &ta, Rotation3::identity(), epsilon = 1.0e-7) &&
+            relative_eq!(  a *  ia, Rotation3::identity(), epsilon = 1.0e-7)
         }
-    }
 
-    fn rotation_between_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
-        if !relative_eq!(a.angle(&b), 0.0, epsilon = 1.0e-7) {
-            let r = Rotation3::rotation_between(&a, &b).unwrap();
-            relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7)
+        fn rotation_inv_2(a: Rotation2<f64>) -> bool {
+            let ta = a.transpose();
+            let ia = a.inverse();
+
+            ta == ia &&
+            relative_eq!(&ta * &a,  Rotation2::identity(), epsilon = 1.0e-7) &&
+            relative_eq!(&ia *  a,  Rotation2::identity(), epsilon = 1.0e-7) &&
+            relative_eq!(  a * &ta, Rotation2::identity(), epsilon = 1.0e-7) &&
+            relative_eq!(  a *  ia, Rotation2::identity(), epsilon = 1.0e-7)
         }
-        else {
-            true
+
+        /*
+         *
+         * Angle between vectors.
+         *
+         */
+        fn angle_is_commutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
+            a.angle(&b) == b.angle(&a)
+        }
+
+        fn angle_is_commutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
+            a.angle(&b) == b.angle(&a)
+        }
+
+        /*
+         *
+         * Rotation matrix between vectors.
+         *
+         */
+        fn rotation_between_is_anticommutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
+            let rab = Rotation2::rotation_between(&a, &b);
+            let rba = Rotation2::rotation_between(&b, &a);
+
+            relative_eq!(rab * rba, Rotation2::identity())
+        }
+
+        fn rotation_between_is_anticommutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
+            let rots = (Rotation3::rotation_between(&a, &b), Rotation3::rotation_between(&b, &a));
+            if let (Some(rab), Some(rba)) = rots {
+                relative_eq!(rab * rba, Rotation3::identity(), epsilon = 1.0e-7)
+            }
+            else {
+                true
+            }
+        }
+
+        fn rotation_between_is_identity(v2: Vector2<f64>, v3: Vector3<f64>) -> bool {
+            let vv2 = 3.42 * v2;
+            let vv3 = 4.23 * v3;
+
+            relative_eq!(v2.angle(&vv2), 0.0, epsilon = 1.0e-7) &&
+            relative_eq!(v3.angle(&vv3), 0.0, epsilon = 1.0e-7) &&
+            relative_eq!(Rotation2::rotation_between(&v2, &vv2), Rotation2::identity()) &&
+            Rotation3::rotation_between(&v3, &vv3).unwrap() == Rotation3::identity()
+        }
+
+        fn rotation_between_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
+            if !relative_eq!(a.angle(&b), 0.0, epsilon = 1.0e-7) {
+                let r = Rotation2::rotation_between(&a, &b);
+                relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7)
+            }
+            else {
+                true
+            }
+        }
+
+        fn rotation_between_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
+            if !relative_eq!(a.angle(&b), 0.0, epsilon = 1.0e-7) {
+                let r = Rotation3::rotation_between(&a, &b).unwrap();
+                relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7)
+            }
+            else {
+                true
+            }
         }
-    }
 
 
-    /*
-     *
-     * Rotation construction.
-     *
-     */
-    fn new_rotation_2(angle: f64) -> bool {
-        let r = Rotation2::new(angle);
+        /*
+         *
+         * Rotation construction.
+         *
+         */
+        fn new_rotation_2(angle: f64) -> bool {
+            let r = Rotation2::new(angle);
 
-        let angle = na::wrap(angle, -f64::pi(), f64::pi());
-        relative_eq!(r.angle(), angle, epsilon = 1.0e-7)
-    }
-
-    fn new_rotation_3(axisangle: Vector3<f64>) -> bool {
-        let r = Rotation3::new(axisangle);
-
-        if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, 0.0) {
             let angle = na::wrap(angle, -f64::pi(), f64::pi());
-            (relative_eq!(r.angle(), angle, epsilon = 1.0e-7) &&
-             relative_eq!(r.axis().unwrap(), axis, epsilon = 1.0e-7)) ||
-            (relative_eq!(r.angle(), -angle, epsilon = 1.0e-7) &&
-             relative_eq!(r.axis().unwrap(), -axis, epsilon = 1.0e-7))
+            relative_eq!(r.angle(), angle, epsilon = 1.0e-7)
         }
-        else {
-            r == Rotation3::identity()
+
+        fn new_rotation_3(axisangle: Vector3<f64>) -> bool {
+            let r = Rotation3::new(axisangle);
+
+            if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, 0.0) {
+                let angle = na::wrap(angle, -f64::pi(), f64::pi());
+                (relative_eq!(r.angle(), angle, epsilon = 1.0e-7) &&
+                 relative_eq!(r.axis().unwrap(), axis, epsilon = 1.0e-7)) ||
+                (relative_eq!(r.angle(), -angle, epsilon = 1.0e-7) &&
+                 relative_eq!(r.axis().unwrap(), -axis, epsilon = 1.0e-7))
+            }
+            else {
+                r == Rotation3::identity()
+            }
         }
-    }
 
-    /*
-     *
-     * Rotation pow.
-     *
-     */
-    fn powf_rotation_2(angle: f64, pow: f64) -> bool {
-        let r = Rotation2::new(angle).powf(pow);
+        /*
+         *
+         * Rotation pow.
+         *
+         */
+        fn powf_rotation_2(angle: f64, pow: f64) -> bool {
+            let r = Rotation2::new(angle).powf(pow);
 
-        let angle  = na::wrap(angle, -f64::pi(), f64::pi());
-        let pangle = na::wrap(angle * pow, -f64::pi(), f64::pi());
-        relative_eq!(r.angle(), pangle, epsilon = 1.0e-7)
-    }
-
-    fn powf_rotation_3(axisangle: Vector3<f64>, pow: f64) -> bool {
-        let r = Rotation3::new(axisangle).powf(pow);
-
-        if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, 0.0) {
-            let angle = na::wrap(angle, -f64::pi(), f64::pi());
+            let angle  = na::wrap(angle, -f64::pi(), f64::pi());
             let pangle = na::wrap(angle * pow, -f64::pi(), f64::pi());
-
-            (relative_eq!(r.angle(), pangle, epsilon = 1.0e-7) &&
-             relative_eq!(r.axis().unwrap(), axis, epsilon = 1.0e-7)) ||
-            (relative_eq!(r.angle(), -pangle, epsilon = 1.0e-7) &&
-             relative_eq!(r.axis().unwrap(), -axis, epsilon = 1.0e-7))
+            relative_eq!(r.angle(), pangle, epsilon = 1.0e-7)
         }
-        else {
-            r == Rotation3::identity()
+
+        fn powf_rotation_3(axisangle: Vector3<f64>, pow: f64) -> bool {
+            let r = Rotation3::new(axisangle).powf(pow);
+
+            if let Some((axis, angle)) = Unit::try_new_and_get(axisangle, 0.0) {
+                let angle = na::wrap(angle, -f64::pi(), f64::pi());
+                let pangle = na::wrap(angle * pow, -f64::pi(), f64::pi());
+
+                (relative_eq!(r.angle(), pangle, epsilon = 1.0e-7) &&
+                 relative_eq!(r.axis().unwrap(), axis, epsilon = 1.0e-7)) ||
+                (relative_eq!(r.angle(), -pangle, epsilon = 1.0e-7) &&
+                 relative_eq!(r.axis().unwrap(), -axis, epsilon = 1.0e-7))
+            }
+            else {
+                r == Rotation3::identity()
+            }
         }
     }
-);
+}
diff --git a/tests/geometry/similarity.rs b/tests/geometry/similarity.rs
index 2364e992..feb8ba1b 100644
--- a/tests/geometry/similarity.rs
+++ b/tests/geometry/similarity.rs
@@ -1,9 +1,9 @@
+#![cfg(feature = "arbitrary")]
 #![allow(non_snake_case)]
 
 use alga::linear::{Transformation, ProjectiveTransformation};
 use na::{Vector3, Point3, Similarity3, Translation3, Isometry3, UnitQuaternion};
 
-#[cfg(feature = "arbitrary")]
 quickcheck!(
     fn inverse_is_identity(i: Similarity3<f64>, p: Point3<f64>, v: Vector3<f64>) -> bool {
         let ii = i.inverse();
diff --git a/tests/geometry/unit_complex.rs b/tests/geometry/unit_complex.rs
index 420b03c1..65837d8b 100644
--- a/tests/geometry/unit_complex.rs
+++ b/tests/geometry/unit_complex.rs
@@ -1,9 +1,8 @@
+#![cfg(feature = "arbitrary")]
 #![allow(non_snake_case)]
 
 use na::{Unit, UnitComplex, Vector2, Point2, Rotation2};
 
-
-#[cfg(feature = "arbitrary")]
 quickcheck!(
 
     /*
diff --git a/tests/linalg/balancing.rs b/tests/linalg/balancing.rs
index 97af9ab6..73fd8230 100644
--- a/tests/linalg/balancing.rs
+++ b/tests/linalg/balancing.rs
@@ -1,9 +1,10 @@
+#![cfg(feature = "arbitrary")]
+
 use std::cmp;
 
 use na::{DMatrix, Matrix4};
 use na::balancing;
 
-#[cfg(feature = "arbitrary")]
 quickcheck! {
     fn balancing_parlett_reinsch(n: usize) -> bool {
         let n = cmp::min(n, 10);
diff --git a/tests/linalg/bidiagonal.rs b/tests/linalg/bidiagonal.rs
index 30eddd4c..a96ce47c 100644
--- a/tests/linalg/bidiagonal.rs
+++ b/tests/linalg/bidiagonal.rs
@@ -1,7 +1,7 @@
+#![cfg(feature = "arbitrary")]
+
 use na::{DMatrix, Matrix2, Matrix4, Matrix5x3, Matrix3x5};
 
-
-#[cfg(feature = "arbitrary")]
 quickcheck! {
     fn bidiagonal(m: DMatrix<f64>) -> bool {
         if m.len() == 0  {
diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs
index 27c9b583..ceb44bdf 100644
--- a/tests/linalg/eigen.rs
+++ b/tests/linalg/eigen.rs
@@ -1,58 +1,61 @@
-use std::cmp;
-
-use na::{DMatrix, Matrix2, Matrix3, Matrix4};
+use na::DMatrix;
 
 #[cfg(feature = "arbitrary")]
-quickcheck! {
-    fn symmetric_eigen(n: usize) -> bool {
-        let n      = cmp::max(1, cmp::min(n, 10));
-        let m      = DMatrix::<f64>::new_random(n, n);
-        let eig    = m.clone().symmetric_eigen();
-        let recomp = eig.recompose();
+mod quickcheck_tests {
+    use std::cmp;
+    use na::{DMatrix, Matrix2, Matrix3, Matrix4};
 
-        println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
+    quickcheck! {
+        fn symmetric_eigen(n: usize) -> bool {
+            let n      = cmp::max(1, cmp::min(n, 10));
+            let m      = DMatrix::<f64>::new_random(n, n);
+            let eig    = m.clone().symmetric_eigen();
+            let recomp = eig.recompose();
 
-        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
-    }
+            println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
 
-    fn symmetric_eigen_singular(n: usize) -> bool {
-        let n      = cmp::max(1, cmp::min(n, 10));
-        let mut m  = DMatrix::<f64>::new_random(n, n);
-        m.row_mut(n / 2).fill(0.0);
-        m.column_mut(n / 2).fill(0.0);
-        let eig    = m.clone().symmetric_eigen();
-        let recomp = eig.recompose();
+            relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+        }
 
-        println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
+        fn symmetric_eigen_singular(n: usize) -> bool {
+            let n      = cmp::max(1, cmp::min(n, 10));
+            let mut m  = DMatrix::<f64>::new_random(n, n);
+            m.row_mut(n / 2).fill(0.0);
+            m.column_mut(n / 2).fill(0.0);
+            let eig    = m.clone().symmetric_eigen();
+            let recomp = eig.recompose();
 
-        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
-    }
+            println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
 
-    fn symmetric_eigen_static_square_4x4(m: Matrix4<f64>) -> bool {
-        let eig    = m.symmetric_eigen();
-        let recomp = eig.recompose();
+            relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+        }
 
-        println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
+        fn symmetric_eigen_static_square_4x4(m: Matrix4<f64>) -> bool {
+            let eig    = m.symmetric_eigen();
+            let recomp = eig.recompose();
 
-        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
-    }
+            println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
 
-    fn symmetric_eigen_static_square_3x3(m: Matrix3<f64>) -> bool {
-        let eig    = m.symmetric_eigen();
-        let recomp = eig.recompose();
+            relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+        }
 
-        println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
+        fn symmetric_eigen_static_square_3x3(m: Matrix3<f64>) -> bool {
+            let eig    = m.symmetric_eigen();
+            let recomp = eig.recompose();
 
-        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
-    }
+            println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
 
-    fn symmetric_eigen_static_square_2x2(m: Matrix2<f64>) -> bool {
-        let eig    = m.symmetric_eigen();
-        let recomp = eig.recompose();
+            relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+        }
 
-        println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
+        fn symmetric_eigen_static_square_2x2(m: Matrix2<f64>) -> bool {
+            let eig    = m.symmetric_eigen();
+            let recomp = eig.recompose();
 
-        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+            println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
+
+            relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+        }
     }
 }
 
diff --git a/tests/linalg/full_piv_lu.rs b/tests/linalg/full_piv_lu.rs
index 2cfa2c29..0f4b0a72 100644
--- a/tests/linalg/full_piv_lu.rs
+++ b/tests/linalg/full_piv_lu.rs
@@ -1,7 +1,4 @@
-use std::cmp;
-use na::{DMatrix, Matrix3, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5,
-         DVector, Vector4};
-
+use na::Matrix3;
 
 #[test]
 fn full_piv_lu_simple() {
@@ -42,114 +39,119 @@ fn full_piv_lu_simple_with_pivot() {
 }
 
 #[cfg(feature = "arbitrary")]
-quickcheck! {
-    fn full_piv_lu(m: DMatrix<f64>) -> bool {
-        let mut m = m;
-        if m.len() == 0 {
-            m = DMatrix::new_random(1, 1);
-        }
+mod quickcheck_tests {
+    use std::cmp;
+    use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
 
-        let lu = m.clone().full_piv_lu();
-        let (p, l, u, q) = lu.unpack();
-        let mut lu = l * u;
-        p.inv_permute_rows(&mut lu);
-        q.inv_permute_columns(&mut lu);
-
-        relative_eq!(m, lu, epsilon = 1.0e-7)
-    }
-
-    fn full_piv_lu_static_3_5(m: Matrix3x5<f64>) -> bool {
-        let lu = m.full_piv_lu();
-        let (p, l, u, q) = lu.unpack();
-        let mut lu = l * u;
-        p.inv_permute_rows(&mut lu);
-        q.inv_permute_columns(&mut lu);
-
-        relative_eq!(m, lu, epsilon = 1.0e-7)
-    }
-
-    fn full_piv_lu_static_5_3(m: Matrix5x3<f64>) -> bool {
-        let lu = m.full_piv_lu();
-        let (p, l, u, q) = lu.unpack();
-        let mut lu = l * u;
-        p.inv_permute_rows(&mut lu);
-        q.inv_permute_columns(&mut lu);
-
-        relative_eq!(m, lu, epsilon = 1.0e-7)
-    }
-
-    fn full_piv_lu_static_square(m: Matrix4<f64>) -> bool {
-        let lu = m.full_piv_lu();
-        let (p, l, u, q) = lu.unpack();
-        let mut lu = l * u;
-        p.inv_permute_rows(&mut lu);
-        q.inv_permute_columns(&mut lu);
-
-        relative_eq!(m, lu, epsilon = 1.0e-7)
-    }
-
-    fn full_piv_lu_solve(n: usize, nb: usize) -> bool {
-        if n != 0 && nb != 0 {
-            let n  = cmp::min(n, 50);  // To avoid slowing down the test too much.
-            let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
-            let m  = DMatrix::<f64>::new_random(n, n);
+    quickcheck! {
+        fn full_piv_lu(m: DMatrix<f64>) -> bool {
+            let mut m = m;
+            if m.len() == 0 {
+                m = DMatrix::new_random(1, 1);
+            }
 
             let lu = m.clone().full_piv_lu();
-            let b1 = DVector::new_random(n);
-            let b2 = DMatrix::new_random(n, nb);
+            let (p, l, u, q) = lu.unpack();
+            let mut lu = l * u;
+            p.inv_permute_rows(&mut lu);
+            q.inv_permute_columns(&mut lu);
 
-            let sol1 = lu.solve(&b1);
-            let sol2 = lu.solve(&b2);
-
-            return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
-                   (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
+            relative_eq!(m, lu, epsilon = 1.0e-7)
         }
 
-        return true;
-    }
+        fn full_piv_lu_static_3_5(m: Matrix3x5<f64>) -> bool {
+            let lu = m.full_piv_lu();
+            let (p, l, u, q) = lu.unpack();
+            let mut lu = l * u;
+            p.inv_permute_rows(&mut lu);
+            q.inv_permute_columns(&mut lu);
 
-    fn full_piv_lu_solve_static(m: Matrix4<f64>) -> bool {
-         let lu = m.full_piv_lu();
-         let b1 = Vector4::new_random();
-         let b2 = Matrix4x3::new_random();
+            relative_eq!(m, lu, epsilon = 1.0e-7)
+        }
 
-         let sol1 = lu.solve(&b1);
-         let sol2 = lu.solve(&b2);
+        fn full_piv_lu_static_5_3(m: Matrix5x3<f64>) -> bool {
+            let lu = m.full_piv_lu();
+            let (p, l, u, q) = lu.unpack();
+            let mut lu = l * u;
+            p.inv_permute_rows(&mut lu);
+            q.inv_permute_columns(&mut lu);
 
-         return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
-                (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
-    }
+            relative_eq!(m, lu, epsilon = 1.0e-7)
+        }
 
-    fn full_piv_lu_inverse(n: usize) -> bool {
-        let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
-        let m = DMatrix::<f64>::new_random(n, n);
+        fn full_piv_lu_static_square(m: Matrix4<f64>) -> bool {
+            let lu = m.full_piv_lu();
+            let (p, l, u, q) = lu.unpack();
+            let mut lu = l * u;
+            p.inv_permute_rows(&mut lu);
+            q.inv_permute_columns(&mut lu);
 
-        let mut l = m.lower_triangle();
-        let mut u = m.upper_triangle();
+            relative_eq!(m, lu, epsilon = 1.0e-7)
+        }
 
-        // Ensure the matrix is well conditioned for inversion.
-        l.fill_diagonal(1.0);
-        u.fill_diagonal(1.0);
-        let m = l * u;
+        fn full_piv_lu_solve(n: usize, nb: usize) -> bool {
+            if n != 0 && nb != 0 {
+                let n  = cmp::min(n, 50);  // To avoid slowing down the test too much.
+                let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
+                let m  = DMatrix::<f64>::new_random(n, n);
 
-        let m1  = m.clone().full_piv_lu().try_inverse().unwrap();
-        let id1 = &m  * &m1;
-        let id2 = &m1 * &m;
+                let lu = m.clone().full_piv_lu();
+                let b1 = DVector::new_random(n);
+                let b2 = DMatrix::new_random(n, nb);
 
-        return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
-    }
+                let sol1 = lu.solve(&b1);
+                let sol2 = lu.solve(&b2);
 
-    fn full_piv_lu_inverse_static(m: Matrix4<f64>) -> bool {
-        let lu = m.full_piv_lu();
+                return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
+                       (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
+            }
 
-        if let Some(m1)  = lu.try_inverse() {
+            return true;
+        }
+
+        fn full_piv_lu_solve_static(m: Matrix4<f64>) -> bool {
+             let lu = m.full_piv_lu();
+             let b1 = Vector4::new_random();
+             let b2 = Matrix4x3::new_random();
+
+             let sol1 = lu.solve(&b1);
+             let sol2 = lu.solve(&b2);
+
+             return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
+                    (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
+        }
+
+        fn full_piv_lu_inverse(n: usize) -> bool {
+            let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
+            let m = DMatrix::<f64>::new_random(n, n);
+
+            let mut l = m.lower_triangle();
+            let mut u = m.upper_triangle();
+
+            // Ensure the matrix is well conditioned for inversion.
+            l.fill_diagonal(1.0);
+            u.fill_diagonal(1.0);
+            let m = l * u;
+
+            let m1  = m.clone().full_piv_lu().try_inverse().unwrap();
             let id1 = &m  * &m1;
             let id2 = &m1 * &m;
 
-            id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
+            return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
         }
-        else {
-            true
+
+        fn full_piv_lu_inverse_static(m: Matrix4<f64>) -> bool {
+            let lu = m.full_piv_lu();
+
+            if let Some(m1)  = lu.try_inverse() {
+                let id1 = &m  * &m1;
+                let id2 = &m1 * &m;
+
+                id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
+            }
+            else {
+                true
+            }
         }
     }
 }
diff --git a/tests/linalg/hessenberg.rs b/tests/linalg/hessenberg.rs
index 6fce8579..67206487 100644
--- a/tests/linalg/hessenberg.rs
+++ b/tests/linalg/hessenberg.rs
@@ -1,7 +1,8 @@
+#![cfg(feature = "arbitrary")]
+
 use std::cmp;
 use na::{DMatrix, Matrix2, Matrix4};
 
-
 #[test]
 fn hessenberg_simple() {
     let m = Matrix2::new(1.0, 0.0,
@@ -11,8 +12,6 @@ fn hessenberg_simple() {
     assert!(relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7))
 }
 
-
-#[cfg(feature = "arbitrary")]
 quickcheck! {
     fn hessenberg(n: usize) -> bool {
         let n = cmp::max(1, cmp::min(n, 50));
diff --git a/tests/linalg/lu.rs b/tests/linalg/lu.rs
index 2e8cc0ce..7b823968 100644
--- a/tests/linalg/lu.rs
+++ b/tests/linalg/lu.rs
@@ -1,7 +1,4 @@
-use std::cmp;
-use na::{DMatrix, Matrix3, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5,
-         DVector, Vector4};
-
+use na::Matrix3;
 
 #[test]
 fn lu_simple() {
@@ -40,110 +37,115 @@ fn lu_simple_with_pivot() {
 }
 
 #[cfg(feature = "arbitrary")]
-quickcheck! {
-    fn lu(m: DMatrix<f64>) -> bool {
-        let mut m = m;
-        if m.len() == 0 {
-            m = DMatrix::new_random(1, 1);
-        }
+mod quickcheck_tests {
+    use std::cmp;
+    use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
 
-        let lu = m.clone().lu();
-        let (p, l, u) = lu.unpack();
-        let mut lu = l * u;
-        p.inv_permute_rows(&mut lu);
-
-        relative_eq!(m, lu, epsilon = 1.0e-7)
-    }
-
-    fn lu_static_3_5(m: Matrix3x5<f64>) -> bool {
-        let lu = m.lu();
-        let (p, l, u) = lu.unpack();
-        let mut lu = l * u;
-        p.inv_permute_rows(&mut lu);
-
-        relative_eq!(m, lu, epsilon = 1.0e-7)
-    }
-
-    fn lu_static_5_3(m: Matrix5x3<f64>) -> bool {
-        let lu = m.lu();
-        let (p, l, u) = lu.unpack();
-        let mut lu = l * u;
-        p.inv_permute_rows(&mut lu);
-
-        relative_eq!(m, lu, epsilon = 1.0e-7)
-    }
-
-    fn lu_static_square(m: Matrix4<f64>) -> bool {
-        let lu = m.lu();
-        let (p, l, u) = lu.unpack();
-        let mut lu = l * u;
-        p.inv_permute_rows(&mut lu);
-
-        relative_eq!(m, lu, epsilon = 1.0e-7)
-    }
-
-    fn lu_solve(n: usize, nb: usize) -> bool {
-        if n != 0 && nb != 0 {
-            let n  = cmp::min(n, 50);  // To avoid slowing down the test too much.
-            let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
-            let m  = DMatrix::<f64>::new_random(n, n);
+    quickcheck! {
+        fn lu(m: DMatrix<f64>) -> bool {
+            let mut m = m;
+            if m.len() == 0 {
+                m = DMatrix::new_random(1, 1);
+            }
 
             let lu = m.clone().lu();
-            let b1 = DVector::new_random(n);
-            let b2 = DMatrix::new_random(n, nb);
+            let (p, l, u) = lu.unpack();
+            let mut lu = l * u;
+            p.inv_permute_rows(&mut lu);
 
-            let sol1 = lu.solve(&b1);
-            let sol2 = lu.solve(&b2);
-
-            return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
-                   (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
+            relative_eq!(m, lu, epsilon = 1.0e-7)
         }
 
-        return true;
-    }
+        fn lu_static_3_5(m: Matrix3x5<f64>) -> bool {
+            let lu = m.lu();
+            let (p, l, u) = lu.unpack();
+            let mut lu = l * u;
+            p.inv_permute_rows(&mut lu);
 
-    fn lu_solve_static(m: Matrix4<f64>) -> bool {
-         let lu = m.lu();
-         let b1 = Vector4::new_random();
-         let b2 = Matrix4x3::new_random();
+            relative_eq!(m, lu, epsilon = 1.0e-7)
+        }
 
-         let sol1 = lu.solve(&b1);
-         let sol2 = lu.solve(&b2);
+        fn lu_static_5_3(m: Matrix5x3<f64>) -> bool {
+            let lu = m.lu();
+            let (p, l, u) = lu.unpack();
+            let mut lu = l * u;
+            p.inv_permute_rows(&mut lu);
 
-         return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
-                (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
-    }
+            relative_eq!(m, lu, epsilon = 1.0e-7)
+        }
 
-    fn lu_inverse(n: usize) -> bool {
-        let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
-        let m = DMatrix::<f64>::new_random(n, n);
+        fn lu_static_square(m: Matrix4<f64>) -> bool {
+            let lu = m.lu();
+            let (p, l, u) = lu.unpack();
+            let mut lu = l * u;
+            p.inv_permute_rows(&mut lu);
 
-        let mut l = m.lower_triangle();
-        let mut u = m.upper_triangle();
+            relative_eq!(m, lu, epsilon = 1.0e-7)
+        }
 
-        // Ensure the matrix is well conditioned for inversion.
-        l.fill_diagonal(1.0);
-        u.fill_diagonal(1.0);
-        let m = l * u;
+        fn lu_solve(n: usize, nb: usize) -> bool {
+            if n != 0 && nb != 0 {
+                let n  = cmp::min(n, 50);  // To avoid slowing down the test too much.
+                let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
+                let m  = DMatrix::<f64>::new_random(n, n);
 
-        let m1  = m.clone().lu().try_inverse().unwrap();
-        let id1 = &m  * &m1;
-        let id2 = &m1 * &m;
+                let lu = m.clone().lu();
+                let b1 = DVector::new_random(n);
+                let b2 = DMatrix::new_random(n, nb);
 
-        return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
-    }
+                let sol1 = lu.solve(&b1);
+                let sol2 = lu.solve(&b2);
 
-    fn lu_inverse_static(m: Matrix4<f64>) -> bool {
-        let lu  = m.lu();
+                return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
+                       (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
+            }
 
-        if let Some(m1) = lu.try_inverse() {
+            return true;
+        }
+
+        fn lu_solve_static(m: Matrix4<f64>) -> bool {
+             let lu = m.lu();
+             let b1 = Vector4::new_random();
+             let b2 = Matrix4x3::new_random();
+
+             let sol1 = lu.solve(&b1);
+             let sol2 = lu.solve(&b2);
+
+             return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
+                    (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
+        }
+
+        fn lu_inverse(n: usize) -> bool {
+            let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
+            let m = DMatrix::<f64>::new_random(n, n);
+
+            let mut l = m.lower_triangle();
+            let mut u = m.upper_triangle();
+
+            // Ensure the matrix is well conditioned for inversion.
+            l.fill_diagonal(1.0);
+            u.fill_diagonal(1.0);
+            let m = l * u;
+
+            let m1  = m.clone().lu().try_inverse().unwrap();
             let id1 = &m  * &m1;
             let id2 = &m1 * &m;
 
-            id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
+            return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
         }
-        else {
-            true
+
+        fn lu_inverse_static(m: Matrix4<f64>) -> bool {
+            let lu  = m.lu();
+
+            if let Some(m1) = lu.try_inverse() {
+                let id1 = &m  * &m1;
+                let id2 = &m1 * &m;
+
+                id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
+            }
+            else {
+                true
+            }
         }
     }
 }
diff --git a/tests/linalg/qr.rs b/tests/linalg/qr.rs
index 8bdc3ce8..394a8b42 100644
--- a/tests/linalg/qr.rs
+++ b/tests/linalg/qr.rs
@@ -1,8 +1,9 @@
+#![cfg(feature = "arbitrary")]
+
 use std::cmp;
 use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5,
          DVector, Vector4};
 
-#[cfg(feature = "arbitrary")]
 quickcheck! {
     fn qr(m: DMatrix<f64>) -> bool {
         let qr = m.clone().qr();
diff --git a/tests/linalg/real_schur.rs b/tests/linalg/real_schur.rs
index 3ef32297..40f8a011 100644
--- a/tests/linalg/real_schur.rs
+++ b/tests/linalg/real_schur.rs
@@ -1,5 +1,4 @@
-use std::cmp;
-use na::{DMatrix, Matrix2, Matrix3, Matrix4};
+use na::{DMatrix, Matrix3, Matrix4};
 
 
 #[test]
@@ -15,49 +14,54 @@ fn schur_simpl_mat3() {
 }
 
 #[cfg(feature = "arbitrary")]
-quickcheck! {
-    fn schur(n: usize) -> bool {
-        let n = cmp::max(1, cmp::min(n, 10));
-        let m = DMatrix::<f64>::new_random(n, n);
+mod quickcheck_tests {
+    use std::cmp;
+    use na::{DMatrix, Matrix2, Matrix3, Matrix4};
 
-        let (vecs, vals) = m.clone().real_schur().unpack();
+    quickcheck! {
+        fn schur(n: usize) -> bool {
+            let n = cmp::max(1, cmp::min(n, 10));
+            let m = DMatrix::<f64>::new_random(n, n);
 
-        if !relative_eq!(&vecs * &vals * vecs.transpose(), m, epsilon = 1.0e-7) {
-            println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
+            let (vecs, vals) = m.clone().real_schur().unpack();
+
+            if !relative_eq!(&vecs * &vals * vecs.transpose(), m, epsilon = 1.0e-7) {
+                println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
+            }
+
+            relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)
         }
 
-        relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)
-    }
+        fn schur_static_mat2(m: Matrix2<f64>) -> bool {
+            let (vecs, vals) = m.clone().real_schur().unpack();
 
-    fn schur_static_mat2(m: Matrix2<f64>) -> bool {
-        let (vecs, vals) = m.clone().real_schur().unpack();
-
-        let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
-        if !ok {
-            println!("{:.5}{:.5}", vecs, vals);
-            println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.transpose());
+            let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
+            if !ok {
+                println!("{:.5}{:.5}", vecs, vals);
+                println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.transpose());
+            }
+            ok
         }
-        ok
-    }
 
-    fn schur_static_mat3(m: Matrix3<f64>) -> bool {
-        let (vecs, vals) = m.clone().real_schur().unpack();
+        fn schur_static_mat3(m: Matrix3<f64>) -> bool {
+            let (vecs, vals) = m.clone().real_schur().unpack();
 
-        let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
-        if !ok {
-            println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
+            let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
+            if !ok {
+                println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
+            }
+            ok
         }
-        ok
-    }
 
-    fn schur_static_mat4(m: Matrix4<f64>) -> bool {
-        let (vecs, vals) = m.clone().real_schur().unpack();
+        fn schur_static_mat4(m: Matrix4<f64>) -> bool {
+            let (vecs, vals) = m.clone().real_schur().unpack();
 
-        let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
-        if !ok {
-            println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
+            let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
+            if !ok {
+                println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
+            }
+            ok
         }
-        ok
     }
 }
 
diff --git a/tests/linalg/solve.rs b/tests/linalg/solve.rs
index 3940c287..f960fb0a 100644
--- a/tests/linalg/solve.rs
+++ b/tests/linalg/solve.rs
@@ -1,3 +1,5 @@
+#![cfg(feature = "arbitrary")]
+
 use na::{Matrix4, Matrix4x5};
 
 fn unzero_diagonal(a: &mut Matrix4<f64>) {
@@ -8,7 +10,6 @@ fn unzero_diagonal(a: &mut Matrix4<f64>) {
     }
 }
 
-#[cfg(feature = "arbitrary")]
 quickcheck! {
     fn solve_lower_triangular(a: Matrix4<f64>, b: Matrix4x5<f64>) -> bool {
         let mut a = a;
diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs
index df107694..5b09c7c8 100644
--- a/tests/linalg/svd.rs
+++ b/tests/linalg/svd.rs
@@ -1,142 +1,143 @@
-use std::cmp;
-
-use na::{DMatrix, Matrix2, Matrix3, Matrix4, Matrix6, Matrix5x2, Matrix5x3, Matrix2x5, Matrix3x5,
-         DVector};
-
+use na::{DMatrix, Matrix6};
 
 #[cfg(feature = "arbitrary")]
-quickcheck! {
-    fn svd(m: DMatrix<f64>) -> bool {
-        if m.len() > 0 {
-            let svd = m.clone().svd(true, true);
-            let recomp_m = svd.clone().recompose();
-            let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-            let ds = DMatrix::from_diagonal(&s);
+mod quickcheck_tests {
+    use std::cmp;
+    use na::{DMatrix, Matrix2, Matrix3, Matrix4, Matrix5x2, Matrix5x3, Matrix2x5, Matrix3x5, DVector};
 
-            println!("{}{}", &m, &u * &ds * &v_t);
+    quickcheck! {
+        fn svd(m: DMatrix<f64>) -> bool {
+            if m.len() > 0 {
+                let svd = m.clone().svd(true, true);
+                let recomp_m = svd.clone().recompose();
+                let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
+                let ds = DMatrix::from_diagonal(&s);
 
-            s.iter().all(|e| *e >= 0.0) &&
-            relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) &&
-            relative_eq!(m, recomp_m, epsilon = 1.0e-5)
-        }
-        else {
-            true
-        }
-    }
+                println!("{}{}", &m, &u * &ds * &v_t);
 
-    fn svd_static_5_3(m: Matrix5x3<f64>) -> bool {
-        let svd = m.svd(true, true);
-        let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-        let ds = Matrix3::from_diagonal(&s);
-
-        s.iter().all(|e| *e >= 0.0) &&
-        relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
-        u.is_orthogonal(1.0e-5) &&
-        v_t.is_orthogonal(1.0e-5)
-    }
-
-    fn svd_static_5_2(m: Matrix5x2<f64>) -> bool {
-        let svd = m.svd(true, true);
-        let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-        let ds = Matrix2::from_diagonal(&s);
-
-        s.iter().all(|e| *e >= 0.0) &&
-        relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
-        u.is_orthogonal(1.0e-5) &&
-        v_t.is_orthogonal(1.0e-5)
-    }
-
-    fn svd_static_3_5(m: Matrix3x5<f64>) -> bool {
-        let svd = m.svd(true, true);
-        let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-
-        let ds = Matrix3::from_diagonal(&s);
-
-        s.iter().all(|e| *e >= 0.0) &&
-        relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
-    }
-
-    fn svd_static_2_5(m: Matrix2x5<f64>) -> bool {
-        let svd = m.svd(true, true);
-        let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-        let ds = Matrix2::from_diagonal(&s);
-
-        s.iter().all(|e| *e >= 0.0) &&
-        relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
-    }
-
-    fn svd_static_square(m: Matrix4<f64>) -> bool {
-        let svd = m.svd(true, true);
-        let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-        let ds = Matrix4::from_diagonal(&s);
-
-        s.iter().all(|e| *e >= 0.0) &&
-        relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
-        u.is_orthogonal(1.0e-5) &&
-        v_t.is_orthogonal(1.0e-5)
-    }
-
-    fn svd_static_square_2x2(m: Matrix2<f64>) -> bool {
-        let svd = m.svd(true, true);
-        let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
-        let ds = Matrix2::from_diagonal(&s);
-
-        s.iter().all(|e| *e >= 0.0) &&
-        relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
-        u.is_orthogonal(1.0e-5) &&
-        v_t.is_orthogonal(1.0e-5)
-    }
-
-    fn svd_pseudo_inverse(m: DMatrix<f64>) -> bool {
-        if m.len() > 0 {
-            let svd = m.clone().svd(true, true);
-            let pinv = svd.pseudo_inverse(1.0e-10);
-
-            if m.nrows() > m.ncols() {
-                println!("{}", &pinv * &m);
-                (pinv * m).is_identity(1.0e-5)
+                s.iter().all(|e| *e >= 0.0) &&
+                relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) &&
+                relative_eq!(m, recomp_m, epsilon = 1.0e-5)
             }
             else {
-                println!("{}", &m * &pinv);
-                (m * pinv).is_identity(1.0e-5)
+                true
             }
         }
-        else {
+
+        fn svd_static_5_3(m: Matrix5x3<f64>) -> bool {
+            let svd = m.svd(true, true);
+            let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
+            let ds = Matrix3::from_diagonal(&s);
+
+            s.iter().all(|e| *e >= 0.0) &&
+            relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
+            u.is_orthogonal(1.0e-5) &&
+            v_t.is_orthogonal(1.0e-5)
+        }
+
+        fn svd_static_5_2(m: Matrix5x2<f64>) -> bool {
+            let svd = m.svd(true, true);
+            let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
+            let ds = Matrix2::from_diagonal(&s);
+
+            s.iter().all(|e| *e >= 0.0) &&
+            relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
+            u.is_orthogonal(1.0e-5) &&
+            v_t.is_orthogonal(1.0e-5)
+        }
+
+        fn svd_static_3_5(m: Matrix3x5<f64>) -> bool {
+            let svd = m.svd(true, true);
+            let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
+
+            let ds = Matrix3::from_diagonal(&s);
+
+            s.iter().all(|e| *e >= 0.0) &&
+            relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
+        }
+
+        fn svd_static_2_5(m: Matrix2x5<f64>) -> bool {
+            let svd = m.svd(true, true);
+            let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
+            let ds = Matrix2::from_diagonal(&s);
+
+            s.iter().all(|e| *e >= 0.0) &&
+            relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
+        }
+
+        fn svd_static_square(m: Matrix4<f64>) -> bool {
+            let svd = m.svd(true, true);
+            let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
+            let ds = Matrix4::from_diagonal(&s);
+
+            s.iter().all(|e| *e >= 0.0) &&
+            relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
+            u.is_orthogonal(1.0e-5) &&
+            v_t.is_orthogonal(1.0e-5)
+        }
+
+        fn svd_static_square_2x2(m: Matrix2<f64>) -> bool {
+            let svd = m.svd(true, true);
+            let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
+            let ds = Matrix2::from_diagonal(&s);
+
+            s.iter().all(|e| *e >= 0.0) &&
+            relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
+            u.is_orthogonal(1.0e-5) &&
+            v_t.is_orthogonal(1.0e-5)
+        }
+
+        fn svd_pseudo_inverse(m: DMatrix<f64>) -> bool {
+            if m.len() > 0 {
+                let svd = m.clone().svd(true, true);
+                let pinv = svd.pseudo_inverse(1.0e-10);
+
+                if m.nrows() > m.ncols() {
+                    println!("{}", &pinv * &m);
+                    (pinv * m).is_identity(1.0e-5)
+                }
+                else {
+                    println!("{}", &m * &pinv);
+                    (m * pinv).is_identity(1.0e-5)
+                }
+            }
+            else {
+                true
+            }
+        }
+
+        fn svd_solve(n: usize, nb: usize) -> bool {
+            let n = cmp::max(1, cmp::min(n, 10));
+            let nb = cmp::min(nb, 10);
+            let m  = DMatrix::<f64>::new_random(n, n);
+
+            let svd = m.clone().svd(true, true);
+
+            if svd.rank(1.0e-7) == n {
+                let b1 = DVector::new_random(n);
+                let b2 = DMatrix::new_random(n, nb);
+
+                let sol1 = svd.solve(&b1, 1.0e-7);
+                let sol2 = svd.solve(&b2, 1.0e-7);
+
+                let recomp = svd.recompose();
+                if !relative_eq!(m, recomp, epsilon = 1.0e-6) {
+                    println!("{}{}", m, recomp);
+                }
+
+                if !relative_eq!(&m * &sol1, b1, epsilon = 1.0e-6) {
+                    println!("Problem 1: {:.6}{:.6}", b1, &m * sol1);
+                    return false;
+                }
+                if !relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6) {
+                    println!("Problem 2: {:.6}{:.6}", b2, &m * sol2);
+                    return false;
+                }
+            }
+
             true
         }
     }
-
-    fn svd_solve(n: usize, nb: usize) -> bool {
-        let n = cmp::max(1, cmp::min(n, 10));
-        let nb = cmp::min(nb, 10);
-        let m  = DMatrix::<f64>::new_random(n, n);
-
-        let svd = m.clone().svd(true, true);
-
-        if svd.rank(1.0e-7) == n {
-            let b1 = DVector::new_random(n);
-            let b2 = DMatrix::new_random(n, nb);
-
-            let sol1 = svd.solve(&b1, 1.0e-7);
-            let sol2 = svd.solve(&b2, 1.0e-7);
-
-            let recomp = svd.recompose();
-            if !relative_eq!(m, recomp, epsilon = 1.0e-6) {
-                println!("{}{}", m, recomp);
-            }
-
-            if !relative_eq!(&m * &sol1, b1, epsilon = 1.0e-6) {
-                println!("Problem 1: {:.6}{:.6}", b1, &m * sol1);
-                return false;
-            }
-            if !relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6) {
-                println!("Problem 2: {:.6}{:.6}", b2, &m * sol2);
-                return false;
-            }
-        }
-
-        true
-    }
 }
 
 // Test proposed on the issue #176 of rulinalg.
diff --git a/tests/linalg/tridiagonal.rs b/tests/linalg/tridiagonal.rs
index 2128d983..6db86aef 100644
--- a/tests/linalg/tridiagonal.rs
+++ b/tests/linalg/tridiagonal.rs
@@ -1,9 +1,9 @@
+#![cfg(feature = "arbitrary")]
+
 use std::cmp;
 
 use na::{DMatrix, Matrix2, Matrix4};
 
-
-#[cfg(feature = "arbitrary")]
 quickcheck! {
     fn symm_tridiagonal(n: usize) -> bool {
         let n = cmp::max(1, cmp::min(n, 50));