Use proptest for all nalgebra tests.

.circleci/config.yml

@@ -60,7 +60,7 @@ jobs:
           command: cargo test --features arbitrary --features serde-serialize --features abomonation-serialize --features sparse --features debug --features io --features compare --features libm --features proptest-support --features slow-tests
       - run:
           name: test nalgebra-glm
-          command: cargo test -p nalgebra-glm  --features arbitrary --features serde-serialize --features abomonation-serialize --features sparse --features debug --features io --features compare --features libm --features slow-tests
+          command: cargo test -p nalgebra-glm  --features arbitrary --features serde-serialize --features abomonation-serialize --features sparse --features debug --features io --features compare --features libm --features proptest-support --features slow-tests
       - run:
           name: test nalgebra-sparse
           # Manifest-path is necessary because cargo otherwise won't correctly forward features

Cargo.toml

@@ -23,8 +23,7 @@ path = "src/lib.rs"
 
 [features]
 default         = [ "std" ]
-std             = [ "matrixmultiply", "rand/std", "rand_distr", "simba/std" ]
-stdweb          = [ "rand/stdweb" ]
+std             = [ "matrixmultiply", "rand/std", "rand/std_rng", "rand_distr", "simba/std" ]
 arbitrary       = [ "quickcheck" ]
 serde-serialize = [ "serde", "num-complex/serde" ]
 abomonation-serialize = [ "abomonation" ]
@@ -44,29 +43,29 @@ slow-tests = []
 [dependencies]
 typenum        = "1.12"
 generic-array  = "0.14"
-rand           = { version = "0.7", default-features = false }
+rand           = { version = "0.8", default-features = false }
 num-traits     = { version = "0.2", default-features = false }
 num-complex    = { version = "0.3", default-features = false }
 num-rational   = { version = "0.3", default-features = false }
 approx         = { version = "0.4", default-features = false }
-simba          = { version = "0.3", default-features = false }
+simba          = { version = "0.4", default-features = false }
 alga           = { version = "0.9", default-features = false, optional = true }
-rand_distr     = { version = "0.3", optional = true }
-matrixmultiply = { version = "0.2", optional = true }
+rand_distr     = { version = "0.4", optional = true }
+matrixmultiply = { version = "0.3", optional = true }
 serde          = { version = "1.0", default-features = false, features = [ "derive" ], optional = true }
 abomonation    = { version = "0.7", optional = true }
 mint           = { version = "0.5", optional = true }
-quickcheck     = { version = "0.9", optional = true }
+quickcheck     = { version = "1", optional = true }
 pest           = { version = "2", optional = true }
 pest_derive    = { version = "2", optional = true }
 bytemuck       = { version = "1.5", optional = true }
 matrixcompare-core = { version = "0.1", optional = true }
-proptest = { version = "0.10", optional = true, default-features = false, features = ["std"] }
+proptest = { version = "1", optional = true, default-features = false, features = ["std"] }
 
 [dev-dependencies]
 serde_json = "1.0"
-rand_xorshift = "0.2"
-rand_isaac = "0.2"
+rand_xorshift = "0.3"
+rand_isaac = "0.3"
 ### Uncomment this line before running benchmarks.
 ### We can't just let this uncommented because that would break
 ### compilation for #[no-std] because of the terrible Cargo bug
@@ -75,7 +74,7 @@ rand_isaac = "0.2"
 
 # For matrix comparison macro
 matrixcompare = "0.2.0"
-itertools = "0.9"
+itertools = "0.10"
 
 [workspace]
 members = [ "nalgebra-lapack", "nalgebra-glm", "nalgebra-sparse" ]

nalgebra-glm/Cargo.toml

@@ -19,7 +19,6 @@ maintenance = { status = "actively-developed" }
 [features]
 default         = [ "std" ]
 std             = [ "nalgebra/std", "simba/std" ]
-stdweb          = [ "nalgebra/stdweb" ]
 arbitrary       = [ "nalgebra/arbitrary" ]
 serde-serialize = [ "nalgebra/serde-serialize" ]
 abomonation-serialize = [ "nalgebra/abomonation-serialize" ]
@@ -27,5 +26,5 @@ abomonation-serialize = [ "nalgebra/abomonation-serialize" ]
 [dependencies]
 num-traits = { version = "0.2", default-features = false }
 approx = { version = "0.4", default-features = false }
-simba = { version = "0.3", default-features = false }
+simba = { version = "0.4", default-features = false }
 nalgebra = { path =  "..", version = "0.24", default-features = false }

nalgebra-lapack/Cargo.toml

@@ -29,16 +29,16 @@ intel-mkl  = ["lapack-src/intel-mkl"]
 [dependencies]
 nalgebra      = { version = "0.24", path = ".." }
 num-traits    = "0.2"
-num-complex   = { version = "0.2", default-features = false }
-simba         = "0.2"
+num-complex   = { version = "0.3", default-features = false }
+simba         = "0.4"
 serde         = { version = "1.0", optional = true }
 serde_derive  = { version = "1.0", optional = true }
-lapack        = { version = "0.16", default-features = false }
-lapack-src    = { version = "0.5", default-features = false }
+lapack        = { version = "0.17", default-features = false }
+lapack-src    = { version = "0.6", default-features = false }
 # clippy = "*"
 
 [dev-dependencies]
 nalgebra   = { version = "0.24", features = [ "arbitrary" ], path = ".." }
-quickcheck = "0.9"
-approx     = "0.3"
-rand       = "0.7"
+quickcheck = "1"
+approx     = "0.4"
+rand       = "0.8"

nalgebra-sparse/Cargo.toml

@@ -14,11 +14,11 @@ slow-tests = []
 [dependencies]
 nalgebra = { version="0.24", path = "../" }
 num-traits = { version = "0.2", default-features = false }
-proptest = { version = "0.10", optional = true }
+proptest = { version = "1.0", optional = true }
 matrixcompare-core = { version = "0.1.0", optional = true }
 
 [dev-dependencies]
-itertools = "0.9"
+itertools = "0.10"
 matrixcompare = { version = "0.2.0", features = [ "proptest-support" ] }
 nalgebra = { version="0.24", path = "../", features = ["compare"] }
 

src/base/construction.rs

@@ -824,8 +824,8 @@ where
 {
     #[inline]
     fn sample<'a, G: Rng + ?Sized>(&self, rng: &'a mut G) -> MatrixMN<N, R, C> {
-        let nrows = R::try_to_usize().unwrap_or_else(|| rng.gen_range(0, 10));
-        let ncols = C::try_to_usize().unwrap_or_else(|| rng.gen_range(0, 10));
+        let nrows = R::try_to_usize().unwrap_or_else(|| rng.gen_range(0..10));
+        let ncols = C::try_to_usize().unwrap_or_else(|| rng.gen_range(0..10));
 
         MatrixMN::from_fn_generic(R::from_usize(nrows), C::from_usize(ncols), |_, _| rng.gen())
     }
@@ -841,9 +841,9 @@ where
     Owned<N, R, C>: Clone + Send,
 {
     #[inline]
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
-        let nrows = R::try_to_usize().unwrap_or(g.gen_range(0, 10));
-        let ncols = C::try_to_usize().unwrap_or(g.gen_range(0, 10));
+    fn arbitrary(g: &mut Gen) -> Self {
+        let nrows = R::try_to_usize().unwrap_or(usize::arbitrary(g) % 10);
+        let ncols = C::try_to_usize().unwrap_or(usize::arbitrary(g) % 10);
 
         Self::from_fn_generic(R::from_usize(nrows), C::from_usize(ncols), |_, _| {
             N::arbitrary(g)

src/base/helper.rs

@@ -7,7 +7,7 @@ use rand::Rng;
 #[cfg(feature = "arbitrary")]
 #[doc(hidden)]
 #[inline]
-pub fn reject<G: Gen, F: FnMut(&T) -> bool, T: Arbitrary>(g: &mut G, f: F) -> T {
+pub fn reject<F: FnMut(&T) -> bool, T: Arbitrary>(g: &mut Gen, f: F) -> T {
     use std::iter;
     iter::repeat(())
         .map(|_| Arbitrary::arbitrary(g))

src/debug/random_orthogonal.rs

@@ -48,9 +48,8 @@ where
     DefaultAllocator: Allocator<N, D, D>,
     Owned<N, D, D>: Clone + Send,
 {
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
-        use rand::Rng;
-        let dim = D::try_to_usize().unwrap_or(g.gen_range(1, 50));
+    fn arbitrary(g: &mut Gen) -> Self {
+        let dim = D::try_to_usize().unwrap_or(1 + usize::arbitrary(g) % 50);
         Self::new(D::from_usize(dim), || N::arbitrary(g))
     }
 }

src/debug/random_sdp.rs

@@ -51,9 +51,8 @@ where
     DefaultAllocator: Allocator<N, D, D>,
     Owned<N, D, D>: Clone + Send,
 {
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
-        use rand::Rng;
-        let dim = D::try_to_usize().unwrap_or(g.gen_range(1, 50));
+    fn arbitrary(g: &mut Gen) -> Self {
+        let dim = D::try_to_usize().unwrap_or(1 + usize::arbitrary(g) % 50);
         Self::new(D::from_usize(dim), || N::arbitrary(g))
     }
 }

src/geometry/dual_quaternion_construction.rs

@@ -108,7 +108,7 @@ where
     N::Element: SimdRealField,
 {
     #[inline]
-    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+    fn arbitrary(rng: &mut Gen) -> Self {
         Self::from_real_and_dual(Arbitrary::arbitrary(rng), Arbitrary::arbitrary(rng))
     }
 }
@@ -216,7 +216,7 @@ where
     N::Element: SimdRealField,
 {
     #[inline]
-    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+    fn arbitrary(rng: &mut Gen) -> Self {
         Self::new_normalize(Arbitrary::arbitrary(rng))
     }
 }

src/geometry/isometry_construction.rs

@@ -102,7 +102,7 @@ where
     DefaultAllocator: Allocator<N, D>,
 {
     #[inline]
-    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+    fn arbitrary(rng: &mut Gen) -> Self {
         Self::from_parts(Arbitrary::arbitrary(rng), Arbitrary::arbitrary(rng))
     }
 }

src/geometry/orthographic.rs

@@ -705,7 +705,7 @@ impl<N: RealField + Arbitrary> Arbitrary for Orthographic3<N>
 where
     Matrix4<N>: Send,
 {
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         let left = Arbitrary::arbitrary(g);
         let right = helper::reject(g, |x: &N| *x > left);
         let bottom = Arbitrary::arbitrary(g);

src/geometry/perspective.rs

@@ -283,7 +283,7 @@ where
 
 #[cfg(feature = "arbitrary")]
 impl<N: RealField + Arbitrary> Arbitrary for Perspective3<N> {
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         let znear = Arbitrary::arbitrary(g);
         let zfar = helper::reject(g, |&x: &N| !(x - znear).is_zero());
         let aspect = helper::reject(g, |&x: &N| !x.is_zero());

src/geometry/point_construction.rs

@@ -156,7 +156,7 @@ where
     <DefaultAllocator as Allocator<N, D>>::Buffer: Send,
 {
     #[inline]
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         Self::from(VectorN::arbitrary(g))
     }
 }

src/geometry/quaternion_construction.rs

@@ -160,7 +160,7 @@ where
     Owned<N, U4>: Send,
 {
     #[inline]
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         Self::new(
             N::arbitrary(g),
             N::arbitrary(g),
@@ -845,7 +845,7 @@ where
     Owned<N, U3>: Send,
 {
     #[inline]
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         let axisangle = Vector3::arbitrary(g);
         Self::from_scaled_axis(axisangle)
     }

src/geometry/rotation_specialization.rs

@@ -275,7 +275,7 @@ where
     Owned<N, U2, U2>: Send,
 {
     #[inline]
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         Self::new(N::arbitrary(g))
     }
 }
@@ -961,7 +961,7 @@ where
     Owned<N, U3>: Send,
 {
     #[inline]
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         Self::new(VectorN::arbitrary(g))
     }
 }

src/geometry/similarity_construction.rs

@@ -114,7 +114,7 @@ where
     Owned<N, D>: Send,
 {
     #[inline]
-    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+    fn arbitrary(rng: &mut Gen) -> Self {
         let mut s: N = Arbitrary::arbitrary(rng);
         while s.is_zero() {
             s = Arbitrary::arbitrary(rng)

src/geometry/translation_construction.rs

@@ -61,13 +61,13 @@ where
 }
 
 #[cfg(feature = "arbitrary")]
-impl<N: Scalar + Arbitrary, D: DimName> Arbitrary for Translation<N, D>
+impl<N: Scalar + Arbitrary + Send, D: DimName> Arbitrary for Translation<N, D>
 where
     DefaultAllocator: Allocator<N, D>,
     Owned<N, D>: Send,
 {
     #[inline]
-    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+    fn arbitrary(rng: &mut Gen) -> Self {
         let v: VectorN<N, D> = Arbitrary::arbitrary(rng);
         Self::from(v)
     }

src/geometry/unit_complex_construction.rs

@@ -395,7 +395,7 @@ where
     N::Element: SimdRealField,
 {
     #[inline]
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         UnitComplex::from_angle(N::arbitrary(g))
     }
 }

src/linalg/decomposition.rs

@@ -154,7 +154,7 @@ impl<N: ComplexField, D: Dim, S: Storage<N, D, D>> Matrix<N, D, D, S> {
     ///
     /// The input matrix `self` is assumed to be symmetric and this decomposition will only read
     /// the upper-triangular part of `self`.
-    pub fn udu(self) -> UDU<N, D>
+    pub fn udu(self) -> Option<UDU<N, D>>
     where
         N: RealField,
         DefaultAllocator: Allocator<N, D> + Allocator<N, D, D>,

src/linalg/udu.rs

@@ -48,7 +48,7 @@ where
     /// the upper-triangular part of `p`.
     ///
     /// Ref.: "Optimal control and estimation-Dover Publications", Robert F. Stengel, (1994) page 360
-    pub fn new(p: MatrixN<N, D>) -> Self {
+    pub fn new(p: MatrixN<N, D>) -> Option<Self> {
         let n = p.ncols();
         let n_dim = p.data.shape().1;
 
@@ -56,6 +56,11 @@ where
         let mut u = MatrixN::zeros_generic(n_dim, n_dim);
 
         d[n - 1] = p[(n - 1, n - 1)];
+
+        if d[n - 1].is_zero() {
+            return None;
+        }
+
         u.column_mut(n - 1)
             .axpy(N::one() / d[n - 1], &p.column(n - 1), N::zero());
 
@@ -67,6 +72,10 @@ where
 
             d[j] = p[(j, j)] - d_j;
 
+            if d[j].is_zero() {
+                return None;
+            }
+
             for i in (0..=j).rev() {
                 let mut u_ij = u[(i, j)];
                 for k in j + 1..n {
@@ -79,7 +88,7 @@ where
             u[(j, j)] = N::one();
         }
 
-        Self { u, d }
+        Some(Self { u, d })
     }
 
     /// Returns the diagonal elements as a matrix

src/proptest/mod.rs

@@ -408,7 +408,7 @@ where
 }
 
 /// A strategy for generating matrices.
-#[derive(Debug)]
+#[derive(Debug, Clone)]
 pub struct MatrixStrategy<NStrategy, R: Dim, C: Dim>
 where
     NStrategy: Strategy,

tests/core/blas.rs

@@ -21,19 +21,20 @@ fn gemm_noncommutative() {
     assert_eq!(res, Matrix2::zero());
 }
 
-#[cfg(feature = "arbitrary")]
-mod blas_quickcheck {
+#[cfg(feature = "proptest-support")]
+mod blas_proptest {
+    use crate::proptest::{PROPTEST_F64, PROPTEST_MATRIX_DIM};
     use na::{DMatrix, DVector};
-    use std::cmp;
+    use proptest::{prop_assert, proptest};
 
-    quickcheck! {
+    proptest! {
         /*
          *
          * Symmetric operators.
          *
          */
-        fn gemv_symm(n: usize, alpha: f64, beta: f64) -> bool {
-            let n = cmp::max(1, cmp::min(n, 50));
+        #[test]
+        fn gemv_symm(n in PROPTEST_MATRIX_DIM, alpha in PROPTEST_F64, beta in PROPTEST_F64) {
             let a = DMatrix::<f64>::new_random(n, n);
             let a = &a * a.transpose();
 
@@ -44,18 +45,16 @@ mod blas_quickcheck {
             y1.gemv(alpha, &a, &x, beta);
             y2.sygemv(alpha, &a.lower_triangle(), &x, beta);
 
-            if !relative_eq!(y1, y2, epsilon = 1.0e-10) {
-                return false;
-            }
+            prop_assert!(relative_eq!(y1, y2, epsilon = 1.0e-10));
 
             y1.gemv(alpha, &a, &x, 0.0);
             y2.sygemv(alpha, &a.lower_triangle(), &x, 0.0);
 
-            relative_eq!(y1, y2, epsilon = 1.0e-10)
+            prop_assert!(relative_eq!(y1, y2, epsilon = 1.0e-10))
         }
 
-        fn gemv_tr(n: usize, alpha: f64, beta: f64) -> bool {
-            let n = cmp::max(1, cmp::min(n, 50));
+        #[test]
+        fn gemv_tr(n in PROPTEST_MATRIX_DIM, alpha in PROPTEST_F64, beta in PROPTEST_F64) {
             let a = DMatrix::<f64>::new_random(n, n);
             let x = DVector::new_random(n);
             let mut y1 = DVector::new_random(n);
@@ -64,18 +63,16 @@ mod blas_quickcheck {
             y1.gemv(alpha, &a, &x, beta);
             y2.gemv_tr(alpha, &a.transpose(), &x, beta);
 
-            if !relative_eq!(y1, y2, epsilon = 1.0e-10) {
-                return false;
-            }
+            prop_assert!(relative_eq!(y1, y2, epsilon = 1.0e-10));
 
             y1.gemv(alpha, &a, &x, 0.0);
             y2.gemv_tr(alpha, &a.transpose(), &x, 0.0);
 
-            relative_eq!(y1, y2, epsilon = 1.0e-10)
+            prop_assert!(relative_eq!(y1, y2, epsilon = 1.0e-10))
         }
 
-        fn ger_symm(n: usize, alpha: f64, beta: f64) -> bool {
-            let n = cmp::max(1, cmp::min(n, 50));
+        #[test]
+        fn ger_symm(n in PROPTEST_MATRIX_DIM, alpha in PROPTEST_F64, beta in PROPTEST_F64) {
             let a = DMatrix::<f64>::new_random(n, n);
             let mut a1 = &a * a.transpose();
             let mut a2 = a1.lower_triangle();
@@ -86,18 +83,16 @@ mod blas_quickcheck {
             a1.ger(alpha, &x, &y, beta);
             a2.syger(alpha, &x, &y, beta);
 
-            if !relative_eq!(a1.lower_triangle(), a2) {
-                return false;
-            }
+            prop_assert!(relative_eq!(a1.lower_triangle(), a2));
 
             a1.ger(alpha, &x, &y, 0.0);
             a2.syger(alpha, &x, &y, 0.0);
 
-            relative_eq!(a1.lower_triangle(), a2)
+            prop_assert!(relative_eq!(a1.lower_triangle(), a2))
         }
 
-        fn quadform(n: usize, alpha: f64, beta: f64) -> bool {
-            let n       = cmp::max(1, cmp::min(n, 50));
+        #[test]
+        fn quadform(n in PROPTEST_MATRIX_DIM, alpha in PROPTEST_F64, beta in PROPTEST_F64) {
             let rhs     = DMatrix::<f64>::new_random(6, n);
             let mid     = DMatrix::<f64>::new_random(6, 6);
             let mut res = DMatrix::new_random(n, n);
@@ -106,13 +101,11 @@ mod blas_quickcheck {
 
             res.quadform(alpha, &mid, &rhs, beta);
 
-            println!("{}{}", res, expected);
-
-            relative_eq!(res, expected, epsilon = 1.0e-7)
+            prop_assert!(relative_eq!(res, expected, epsilon = 1.0e-7))
         }
 
-        fn quadform_tr(n: usize, alpha: f64, beta: f64) -> bool {
-            let n       = cmp::max(1, cmp::min(n, 50));
+        #[test]
+        fn quadform_tr(n in PROPTEST_MATRIX_DIM, alpha in PROPTEST_F64, beta in PROPTEST_F64) {
             let lhs     = DMatrix::<f64>::new_random(6, n);
             let mid     = DMatrix::<f64>::new_random(n, n);
             let mut res = DMatrix::new_random(6, 6);
@@ -121,9 +114,7 @@ mod blas_quickcheck {
 
             res.quadform_tr(alpha, &lhs, &mid , beta);
 
-            println!("{}{}", res, expected);
-
-            relative_eq!(res, expected, epsilon = 1.0e-7)
+            prop_assert!(relative_eq!(res, expected, epsilon = 1.0e-7))
         }
     }
 }

tests/core/conversion.rs

@@ -1,44 +1,49 @@
-#![cfg(all(feature = "arbitrary", feature = "alga"))]
+#![cfg(all(feature = "proptest-support", feature = "alga"))]
 use alga::linear::Transformation;
 use na::{
     self, Affine3, Isometry3, Matrix2, Matrix2x3, Matrix2x4, Matrix2x5, Matrix2x6, Matrix3,
     Matrix3x2, Matrix3x4, Matrix3x5, Matrix3x6, Matrix4, Matrix4x2, Matrix4x3, Matrix4x5,
     Matrix4x6, Matrix5, Matrix5x2, Matrix5x3, Matrix5x4, Matrix5x6, Matrix6, Matrix6x2, Matrix6x3,
-    Matrix6x4, Matrix6x5, Point3, Projective3, Rotation3, RowVector1, RowVector2, RowVector3,
-    RowVector4, RowVector5, RowVector6, Similarity3, Transform3, Translation3, UnitQuaternion,
-    Vector1, Vector2, Vector3, Vector4, Vector5, Vector6,
+    Matrix6x4, Matrix6x5, Projective3, Rotation3, RowVector1, RowVector2, RowVector3, RowVector4,
+    RowVector5, RowVector6, Similarity3, Transform3, UnitQuaternion, Vector1, Vector2, Vector3,
+    Vector4, Vector5, Vector6,
 };
 use na::{DMatrix, DMatrixSlice, DMatrixSliceMut, MatrixSlice, MatrixSliceMut};
 use na::{U1, U3, U4};
 
-quickcheck! {
-    fn translation_conversion(t: Translation3<f64>, v: Vector3<f64>, p: Point3<f64>) -> bool {
+use crate::proptest::*;
+use proptest::{prop_assert, prop_assert_eq, proptest};
+
+proptest! {
+    #[test]
+    fn translation_conversion(t in translation3(), v in vector3(), p in point3()) {
         let iso: Isometry3<f64>   = na::convert(t);
         let sim: Similarity3<f64> = na::convert(t);
         let aff: Affine3<f64>     = na::convert(t);
         let prj: Projective3<f64> = na::convert(t);
         let tr:  Transform3<f64>  = na::convert(t);
 
-        t == na::try_convert(iso).unwrap() &&
-        t == na::try_convert(sim).unwrap() &&
-        t == na::try_convert(aff).unwrap() &&
-        t == na::try_convert(prj).unwrap() &&
-        t == na::try_convert(tr).unwrap()  &&
+        prop_assert_eq!(t, na::try_convert(iso).unwrap());
+        prop_assert_eq!(t, na::try_convert(sim).unwrap());
+        prop_assert_eq!(t, na::try_convert(aff).unwrap());
+        prop_assert_eq!(t, na::try_convert(prj).unwrap());
+        prop_assert_eq!(t, na::try_convert(tr).unwrap() );
 
-        t.transform_vector(&v) == iso * v &&
-        t.transform_vector(&v) == sim * v &&
-        t.transform_vector(&v) == aff * v &&
-        t.transform_vector(&v) == prj * v &&
-        t.transform_vector(&v) == tr  * v &&
+        prop_assert_eq!(t.transform_vector(&v), iso * v);
+        prop_assert_eq!(t.transform_vector(&v), sim * v);
+        prop_assert_eq!(t.transform_vector(&v), aff * v);
+        prop_assert_eq!(t.transform_vector(&v), prj * v);
+        prop_assert_eq!(t.transform_vector(&v), tr  * v);
 
-        t * p == iso * p &&
-        t * p == sim * p &&
-        t * p == aff * p &&
-        t * p == prj * p &&
-        t * p == tr  * p
+        prop_assert_eq!(t * p, iso * p);
+        prop_assert_eq!(t * p, sim * p);
+        prop_assert_eq!(t * p, aff * p);
+        prop_assert_eq!(t * p, prj * p);
+        prop_assert_eq!(t * p, tr  * p);
     }
 
-    fn rotation_conversion(r: Rotation3<f64>, v: Vector3<f64>, p: Point3<f64>) -> bool {
+    #[test]
+    fn rotation_conversion(r in rotation3(), v in vector3(), p in point3()) {
         let uq:  UnitQuaternion<f64> = na::convert(r);
         let iso: Isometry3<f64>      = na::convert(r);
         let sim: Similarity3<f64>    = na::convert(r);
@@ -46,30 +51,31 @@ quickcheck! {
         let prj: Projective3<f64>    = na::convert(r);
         let tr:  Transform3<f64>     = na::convert(r);
 
-        relative_eq!(r, na::try_convert(uq).unwrap(),  epsilon = 1.0e-7) &&
-        relative_eq!(r, na::try_convert(iso).unwrap(), epsilon = 1.0e-7) &&
-        relative_eq!(r, na::try_convert(sim).unwrap(), epsilon = 1.0e-7) &&
-        r == na::try_convert(aff).unwrap() &&
-        r == na::try_convert(prj).unwrap() &&
-        r == na::try_convert(tr).unwrap()  &&
+        prop_assert!(relative_eq!(r, na::try_convert(uq).unwrap(),  epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(r, na::try_convert(iso).unwrap(), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(r, na::try_convert(sim).unwrap(), epsilon = 1.0e-7));
+        prop_assert_eq!(r, na::try_convert(aff).unwrap());
+        prop_assert_eq!(r, na::try_convert(prj).unwrap());
+        prop_assert_eq!(r, na::try_convert(tr).unwrap() );
 
         // NOTE: we need relative_eq because Isometry and Similarity use quaternions.
-        relative_eq!(r * v, uq  * v, epsilon = 1.0e-7) &&
-        relative_eq!(r * v, iso * v, epsilon = 1.0e-7) &&
-        relative_eq!(r * v, sim * v, epsilon = 1.0e-7) &&
-        r * v == aff * v &&
-        r * v == prj * v &&
-        r * v == tr  * v &&
+        prop_assert!(relative_eq!(r * v, uq  * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(r * v, iso * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(r * v, sim * v, epsilon = 1.0e-7));
+        prop_assert_eq!(r * v, aff * v);
+        prop_assert_eq!(r * v, prj * v);
+        prop_assert_eq!(r * v, tr  * v);
 
-        relative_eq!(r * p, uq  * p, epsilon = 1.0e-7) &&
-        relative_eq!(r * p, iso * p, epsilon = 1.0e-7) &&
-        relative_eq!(r * p, sim * p, epsilon = 1.0e-7) &&
-        r * p == aff * p &&
-        r * p == prj * p &&
-        r * p == tr  * p
+        prop_assert!(relative_eq!(r * p, uq  * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(r * p, iso * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(r * p, sim * p, epsilon = 1.0e-7));
+        prop_assert_eq!(r * p, aff * p);
+        prop_assert_eq!(r * p, prj * p);
+        prop_assert_eq!(r * p, tr  * p);
     }
 
-    fn unit_quaternion_conversion(uq: UnitQuaternion<f64>, v: Vector3<f64>, p: Point3<f64>) -> bool {
+    #[test]
+    fn unit_quaternion_conversion(uq in unit_quaternion(), v in vector3(), p in point3()) {
         let rot: Rotation3<f64>   = na::convert(uq);
         let iso: Isometry3<f64>   = na::convert(uq);
         let sim: Similarity3<f64> = na::convert(uq);
@@ -77,68 +83,70 @@ quickcheck! {
         let prj: Projective3<f64> = na::convert(uq);
         let tr:  Transform3<f64>  = na::convert(uq);
 
-        uq == na::try_convert(iso).unwrap() &&
-        uq == na::try_convert(sim).unwrap() &&
-        relative_eq!(uq, na::try_convert(rot).unwrap(), epsilon = 1.0e-7) &&
-        relative_eq!(uq, na::try_convert(aff).unwrap(), epsilon = 1.0e-7) &&
-        relative_eq!(uq, na::try_convert(prj).unwrap(), epsilon = 1.0e-7) &&
-        relative_eq!(uq, na::try_convert(tr).unwrap(), epsilon = 1.0e-7)  &&
+        prop_assert_eq!(uq, na::try_convert(iso).unwrap());
+        prop_assert_eq!(uq, na::try_convert(sim).unwrap());
+        prop_assert!(relative_eq!(uq, na::try_convert(rot).unwrap(), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(uq, na::try_convert(aff).unwrap(), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(uq, na::try_convert(prj).unwrap(), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(uq, na::try_convert(tr).unwrap(), epsilon = 1.0e-7) );
 
         // NOTE: iso and sim use unit quaternions for the rotation so conversions to them are exact.
-        relative_eq!(uq * v, rot * v, epsilon = 1.0e-7) &&
-        uq * v == iso * v &&
-        uq * v == sim * v &&
-        relative_eq!(uq * v, aff * v, epsilon = 1.0e-7) &&
-        relative_eq!(uq * v, prj * v, epsilon = 1.0e-7) &&
-        relative_eq!(uq * v, tr  * v, epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!(uq * v, rot * v, epsilon = 1.0e-7));
+        prop_assert_eq!(uq * v, iso * v);
+        prop_assert_eq!(uq * v, sim * v);
+        prop_assert!(relative_eq!(uq * v, aff * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(uq * v, prj * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(uq * v, tr  * v, epsilon = 1.0e-7));
 
-        relative_eq!(uq * p, rot * p, epsilon = 1.0e-7) &&
-        uq * p == iso * p &&
-        uq * p == sim * p &&
-        relative_eq!(uq * p, aff * p, epsilon = 1.0e-7) &&
-        relative_eq!(uq * p, prj * p, epsilon = 1.0e-7) &&
-        relative_eq!(uq * p, tr  * p, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(uq * p, rot * p, epsilon = 1.0e-7));
+        prop_assert_eq!(uq * p, iso * p);
+        prop_assert_eq!(uq * p, sim * p);
+        prop_assert!(relative_eq!(uq * p, aff * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(uq * p, prj * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(uq * p, tr  * p, epsilon = 1.0e-7));
     }
 
-    fn isometry_conversion(iso: Isometry3<f64>, v: Vector3<f64>, p: Point3<f64>) -> bool {
+    #[test]
+    fn isometry_conversion(iso in isometry3(), v in vector3(), p in point3()) {
         let sim: Similarity3<f64> = na::convert(iso);
         let aff: Affine3<f64>     = na::convert(iso);
         let prj: Projective3<f64> = na::convert(iso);
         let tr:  Transform3<f64>  = na::convert(iso);
 
 
-        iso == na::try_convert(sim).unwrap() &&
-        relative_eq!(iso, na::try_convert(aff).unwrap(), epsilon = 1.0e-7) &&
-        relative_eq!(iso, na::try_convert(prj).unwrap(), epsilon = 1.0e-7) &&
-        relative_eq!(iso, na::try_convert(tr).unwrap(), epsilon = 1.0e-7)  &&
+        prop_assert_eq!(iso, na::try_convert(sim).unwrap());
+        prop_assert!(relative_eq!(iso, na::try_convert(aff).unwrap(), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(iso, na::try_convert(prj).unwrap(), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(iso, na::try_convert(tr).unwrap(), epsilon = 1.0e-7) );
 
-        iso * v == sim * v &&
-        relative_eq!(iso * v, aff * v, epsilon = 1.0e-7) &&
-        relative_eq!(iso * v, prj * v, epsilon = 1.0e-7) &&
-        relative_eq!(iso * v, tr  * v, epsilon = 1.0e-7) &&
+        prop_assert_eq!(iso * v, sim * v);
+        prop_assert!(relative_eq!(iso * v, aff * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(iso * v, prj * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(iso * v, tr  * v, epsilon = 1.0e-7));
 
-        iso * p == sim * p &&
-        relative_eq!(iso * p, aff * p, epsilon = 1.0e-7) &&
-        relative_eq!(iso * p, prj * p, epsilon = 1.0e-7) &&
-        relative_eq!(iso * p, tr  * p, epsilon = 1.0e-7)
+        prop_assert_eq!(iso * p, sim * p);
+        prop_assert!(relative_eq!(iso * p, aff * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(iso * p, prj * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(iso * p, tr  * p, epsilon = 1.0e-7));
     }
 
-    fn similarity_conversion(sim: Similarity3<f64>, v: Vector3<f64>, p: Point3<f64>) -> bool {
+    #[test]
+    fn similarity_conversion(sim in similarity3(), v in vector3(), p in point3()) {
         let aff: Affine3<f64>     = na::convert(sim);
         let prj: Projective3<f64> = na::convert(sim);
         let tr:  Transform3<f64>  = na::convert(sim);
 
-        relative_eq!(sim, na::try_convert(aff).unwrap(), epsilon = 1.0e-7) &&
-        relative_eq!(sim, na::try_convert(prj).unwrap(), epsilon = 1.0e-7) &&
-        relative_eq!(sim, na::try_convert(tr).unwrap(),  epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!(sim, na::try_convert(aff).unwrap(), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(sim, na::try_convert(prj).unwrap(), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(sim, na::try_convert(tr).unwrap(),  epsilon = 1.0e-7));
 
-        relative_eq!(sim * v, aff * v, epsilon = 1.0e-7) &&
-        relative_eq!(sim * v, prj * v, epsilon = 1.0e-7) &&
-        relative_eq!(sim * v, tr  * v, epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!(sim * v, aff * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(sim * v, prj * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(sim * v, tr  * v, epsilon = 1.0e-7));
 
-        relative_eq!(sim * p, aff * p, epsilon = 1.0e-7) &&
-        relative_eq!(sim * p, prj * p, epsilon = 1.0e-7) &&
-        relative_eq!(sim * p, tr  * p, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(sim * p, aff * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(sim * p, prj * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(sim * p, tr  * p, epsilon = 1.0e-7));
     }
 
     // XXX test Transform

tests/core/helper.rs

@@ -12,7 +12,7 @@ pub struct RandComplex<N>(pub Complex<N>);
 
 impl<N: Arbitrary + RealField> Arbitrary for RandComplex<N> {
     #[inline]
-    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+    fn arbitrary(rng: &mut Gen) -> Self {
         let im = Arbitrary::arbitrary(rng);
         let re = Arbitrary::arbitrary(rng);
         RandComplex(Complex::new(re, im))
@@ -38,7 +38,7 @@ pub struct RandScalar<N>(pub N);
 
 impl<N: Arbitrary> Arbitrary for RandScalar<N> {
     #[inline]
-    fn arbitrary<G: Gen>(rng: &mut G) -> Self {
+    fn arbitrary(rng: &mut Gen) -> Self {
         RandScalar(Arbitrary::arbitrary(rng))
     }
 }

tests/core/matrix.rs

@@ -829,151 +829,145 @@ fn swizzle() {
     assert_eq!(c.zyz(), Vector3::new(3.0, 2.0, 3.0));
 }
 
-#[cfg(feature = "arbitrary")]
+#[cfg(feature = "proptest-support")]
 mod transposition_tests {
     use super::*;
-    use na::Matrix4x6;
+    use crate::proptest::{dmatrix, matrix, vector4, PROPTEST_F64};
+    use na::{U2, U3, U4, U6};
+    use proptest::{prop_assert, prop_assert_eq, proptest};
 
-    quickcheck! {
-        fn transpose_transpose_is_self(m: Matrix2x3<f64>) -> bool {
-            m.transpose().transpose() == m
+    proptest! {
+        #[test]
+        fn transpose_transpose_is_self(m in matrix(PROPTEST_F64, U2, U3)) {
+            prop_assert_eq!(m.transpose().transpose(), m)
         }
 
-        fn transpose_mut_transpose_mut_is_self(m: Matrix3<f64>) -> bool {
+        #[test]
+        fn transpose_mut_transpose_mut_is_self(m in matrix(PROPTEST_F64, U3, U3)) {
             let mut mm = m;
             mm.transpose_mut();
             mm.transpose_mut();
-            m == mm
+            prop_assert_eq!(m, mm)
         }
 
-        fn transpose_transpose_is_id_dyn(m: DMatrix<f64>) -> bool {
-            m.transpose().transpose() == m
+        #[test]
+        fn transpose_transpose_is_id_dyn(m in dmatrix()) {
+            prop_assert_eq!(m.transpose().transpose(), m)
         }
 
-        fn check_transpose_components_dyn(m: DMatrix<f64>) -> bool {
+        #[test]
+        fn check_transpose_components_dyn(m in dmatrix()) {
             let tr = m.transpose();
             let (nrows, ncols) = m.shape();
 
-            if nrows != tr.shape().1 || ncols != tr.shape().0 {
-                return false
-            }
+            prop_assert!(nrows == tr.shape().1 && ncols == tr.shape().0);
 
             for i in 0 .. nrows {
                 for j in 0 .. ncols {
-                    if m[(i, j)] != tr[(j, i)] {
-                        return false
+                    prop_assert_eq!(m[(i, j)], tr[(j, i)]);
                 }
             }
         }
 
-            true
-        }
-
-        fn tr_mul_is_transpose_then_mul(m: Matrix4x6<f64>, v: Vector4<f64>) -> bool {
-            relative_eq!(m.transpose() * v, m.tr_mul(&v), epsilon = 1.0e-7)
+        #[test]
+        fn tr_mul_is_transpose_then_mul(m in matrix(PROPTEST_F64, U4, U6), v in vector4()) {
+            prop_assert!(relative_eq!(m.transpose() * v, m.tr_mul(&v), epsilon = 1.0e-7))
         }
     }
 }
 
-#[cfg(feature = "arbitrary")]
+#[cfg(feature = "proptest-support")]
 mod inversion_tests {
     use super::*;
+    use crate::proptest::*;
     use na::Matrix1;
+    use proptest::{prop_assert, proptest};
 
-    quickcheck! {
-        fn self_mul_inv_is_id_dim1(m: Matrix1<f64>) -> bool {
+    proptest! {
+        #[test]
+        fn self_mul_inv_is_id_dim1(m in matrix1()) {
             if let Some(im) = m.try_inverse() {
                 let id = Matrix1::one();
-                relative_eq!(im * m, id, epsilon = 1.0e-7) &&
-                relative_eq!(m * im, id, epsilon = 1.0e-7)
-            }
-            else {
-                true
+                prop_assert!(relative_eq!(im * m, id, epsilon = 1.0e-7));
+                prop_assert!(relative_eq!(m * im, id, epsilon = 1.0e-7));
             }
         }
 
-        fn self_mul_inv_is_id_dim2(m: Matrix2<f64>) -> bool {
+        #[test]
+        fn self_mul_inv_is_id_dim2(m in matrix2()) {
             if let Some(im) = m.try_inverse() {
                 let id = Matrix2::one();
-                relative_eq!(im * m, id, epsilon = 1.0e-7) &&
-                relative_eq!(m * im, id, epsilon = 1.0e-7)
-            }
-            else {
-                true
+                prop_assert!(relative_eq!(im * m, id, epsilon = 1.0e-7));
+                prop_assert!(relative_eq!(m * im, id, epsilon = 1.0e-7));
             }
         }
 
-        fn self_mul_inv_is_id_dim3(m: Matrix3<f64>) -> bool {
+        #[test]
+        fn self_mul_inv_is_id_dim3(m in matrix3()) {
             if let Some(im) = m.try_inverse() {
                 let id = Matrix3::one();
-                relative_eq!(im * m, id, epsilon = 1.0e-7) &&
-                relative_eq!(m * im, id, epsilon = 1.0e-7)
-            }
-            else {
-                true
+                prop_assert!(relative_eq!(im * m, id, epsilon = 1.0e-7));
+                prop_assert!(relative_eq!(m * im, id, epsilon = 1.0e-7));
             }
         }
 
-        fn self_mul_inv_is_id_dim4(m: Matrix4<f64>) -> bool {
+        #[test]
+        fn self_mul_inv_is_id_dim4(m in matrix4()) {
             if let Some(im) = m.try_inverse() {
                 let id = Matrix4::one();
-                relative_eq!(im * m, id, epsilon = 1.0e-7) &&
-                relative_eq!(m * im, id, epsilon = 1.0e-7)
-            }
-            else {
-                true
+                prop_assert!(relative_eq!(im * m, id, epsilon = 1.0e-7));
+                prop_assert!(relative_eq!(m * im, id, epsilon = 1.0e-7));
             }
         }
 
-        fn self_mul_inv_is_id_dim6(m: Matrix6<f64>) -> bool {
+        #[test]
+        fn self_mul_inv_is_id_dim6(m in matrix6()) {
             if let Some(im) = m.try_inverse() {
                 let id = Matrix6::one();
-                relative_eq!(im * m, id, epsilon = 1.0e-7) &&
-                relative_eq!(m * im, id, epsilon = 1.0e-7)
-            }
-            else {
-                true
+                prop_assert!(relative_eq!(im * m, id, epsilon = 1.0e-7));
+                prop_assert!(relative_eq!(m * im, id, epsilon = 1.0e-7));
             }
         }
     }
 }
 
-#[cfg(feature = "arbitrary")]
+#[cfg(feature = "proptest-support")]
 mod normalization_tests {
-    use super::*;
+    use crate::proptest::*;
+    use proptest::{prop_assert, proptest};
 
-    quickcheck! {
-        fn normalized_vec_norm_is_one(v: Vector3<f64>) -> bool {
+    proptest! {
+        #[test]
+        fn normalized_vec_norm_is_one(v in vector3()) {
             if let Some(nv) = v.try_normalize(1.0e-10) {
-                relative_eq!(nv.norm(), 1.0, epsilon = 1.0e-7)
-            }
-            else {
-                true
+                prop_assert!(relative_eq!(nv.norm(), 1.0, epsilon = 1.0e-7));
             }
         }
 
-        fn normalized_vec_norm_is_one_dyn(v: DVector<f64>) -> bool {
+        #[test]
+        fn normalized_vec_norm_is_one_dyn(v in dvector()) {
             if let Some(nv) = v.try_normalize(1.0e-10) {
-                relative_eq!(nv.norm(), 1.0, epsilon = 1.0e-7)
-            }
-            else {
-                true
+                prop_assert!(relative_eq!(nv.norm(), 1.0, epsilon = 1.0e-7));
             }
         }
     }
 }
 
-#[cfg(all(feature = "arbitrary", feature = "alga"))]
-// TODO: move this to alga ?
+#[cfg(all(feature = "proptest-support", feature = "alga"))]
+// TODO: move this to alga ?
 mod finite_dim_inner_space_tests {
     use super::*;
+    use crate::proptest::*;
     use alga::linear::FiniteDimInnerSpace;
+    use proptest::collection::vec;
+    use proptest::{prop_assert, proptest};
     use std::fmt::Display;
 
     macro_rules! finite_dim_inner_space_test(
-        ($($Vector: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$(
-            quickcheck!{
-                fn $orthonormal_subspace(vs: Vec<$Vector<f64>>) -> bool {
+        ($($Vector: ident, $vstrategy: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$(
+            proptest! {
+                #[test]
+                fn $orthonormal_subspace(vs in vec($vstrategy(), 0..10)) {
                     let mut given_basis = vs.clone();
                     let given_basis_dim = $Vector::orthonormalize(&mut given_basis[..]);
                     let mut ortho_basis = Vec::new();
@@ -982,29 +976,21 @@ mod finite_dim_inner_space_tests {
                         |e| { ortho_basis.push(*e); true }
                     );
 
-                    if !is_subspace_basis(&ortho_basis[..]) {
-                        return false;
-                    }
+                    prop_assert!(is_subspace_basis(&ortho_basis[..]));
 
                     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;
+                            prop_assert!(relative_eq!(v.dot(b), 0.0, epsilon = 1.0e-7));
                         }
                     }
                 }
 
-                    true
-                }
-
-                fn $orthonormalization(vs: Vec<$Vector<f64>>) -> bool {
+                #[test]
+                fn $orthonormalization(vs in vec($vstrategy(), 0..10)) {
                     let mut basis = vs.clone();
                     let subdim = $Vector::orthonormalize(&mut basis[..]);
 
-                    if !is_subspace_basis(&basis[.. subdim]) {
-                        return false;
-                    }
+                    prop_assert!(is_subspace_basis(&basis[.. subdim]));
 
                     for mut e in vs {
                         for b in &basis[.. subdim] {
@@ -1012,26 +998,20 @@ mod finite_dim_inner_space_tests {
                         }
 
                         // 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;
+                        prop_assert!(relative_eq!(e.norm(), 0.0, epsilon = 1.0e-7));
                     }
                 }
-
-                    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;
+        Vector1, vector1, orthonormal_subspace_basis1, orthonormalize1;
+        Vector2, vector2, orthonormal_subspace_basis2, orthonormalize2;
+        Vector3, vector3, orthonormal_subspace_basis3, orthonormalize3;
+        Vector4, vector4, orthonormal_subspace_basis4, orthonormalize4;
+        Vector5, vector5, orthonormal_subspace_basis5, orthonormalize5;
+        Vector6, vector6, orthonormal_subspace_basis6, orthonormalize6;
     );
 
     /*
@@ -1039,7 +1019,6 @@ mod finite_dim_inner_space_tests {
      * Helper functions.
      *
      */
-    #[cfg(feature = "arbitrary")]
     fn is_subspace_basis<T: FiniteDimInnerSpace<RealField = f64, ComplexField = f64> + Display>(
         vs: &[T],
     ) -> bool {

tests/core/matrixcompare.rs

@@ -4,34 +4,32 @@
 //! The tests here only check that the necessary trait implementations are correctly implemented,
 //! in addition to some sanity checks with example input.
 
+use matrixcompare::assert_matrix_eq;
 use nalgebra::{MatrixMN, U4, U5};
 
-#[cfg(feature = "arbitrary")]
-use nalgebra::DMatrix;
+#[cfg(feature = "proptest-support")]
+use {
+    crate::proptest::*,
+    matrixcompare::DenseAccess,
+    nalgebra::DMatrix,
+    proptest::{prop_assert_eq, proptest},
+};
 
-use matrixcompare::assert_matrix_eq;
-
-#[cfg(feature = "arbitrary")]
-use matrixcompare::DenseAccess;
-
-#[cfg(feature = "arbitrary")]
-quickcheck! {
-    fn fetch_single_is_equivalent_to_index_f64(matrix: DMatrix<f64>) -> bool {
+#[cfg(feature = "proptest-support")]
+proptest! {
+    #[test]
+    fn fetch_single_is_equivalent_to_index_f64(matrix in dmatrix()) {
         for i in 0 .. matrix.nrows() {
             for j in 0 .. matrix.ncols() {
-                if matrix.fetch_single(i, j) != *matrix.index((i, j)) {
-                    return false;
+                prop_assert_eq!(matrix.fetch_single(i, j), *matrix.index((i, j)));
             }
         }
     }
 
-        true
-    }
-
-    fn matrixcompare_shape_agrees_with_matrix(matrix: DMatrix<f64>) -> bool {
-        matrix.nrows() == <DMatrix<f64> as matrixcompare::Matrix<f64>>::rows(&matrix)
-        &&
-        matrix.ncols() == <DMatrix<f64> as matrixcompare::Matrix<f64>>::cols(&matrix)
+    #[test]
+    fn matrixcompare_shape_agrees_with_matrix(matrix in dmatrix()) {
+        prop_assert_eq!(matrix.nrows(), <DMatrix<f64> as matrixcompare::Matrix<f64>>::rows(&matrix));
+        prop_assert_eq!(matrix.ncols(), <DMatrix<f64> as matrixcompare::Matrix<f64>>::cols(&matrix));
     }
 }
 

tests/geometry/dual_quaternion.rs

@@ -1,36 +1,40 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 #![allow(non_snake_case)]
 
-use na::{
-    DualQuaternion, Isometry3, Point3, Translation3, UnitDualQuaternion, UnitQuaternion, Vector3,
-};
+use na::{DualQuaternion, Point3, UnitDualQuaternion, Vector3};
 
-quickcheck!(
-    fn isometry_equivalence(iso: Isometry3<f64>, p: Point3<f64>, v: Vector3<f64>) -> bool {
+use crate::proptest::*;
+use proptest::{prop_assert, proptest};
+
+proptest!(
+    #[test]
+    fn isometry_equivalence(iso in isometry3(), p in point3(), v in vector3()) {
         let dq = UnitDualQuaternion::from_isometry(&iso);
 
-        relative_eq!(iso * p, dq * p, epsilon = 1.0e-7)
-            && relative_eq!(iso * v, dq * v, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(iso * p, dq * p, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(iso * v, dq * v, epsilon = 1.0e-7));
     }
 
-    fn inverse_is_identity(i: UnitDualQuaternion<f64>, p: Point3<f64>, v: Vector3<f64>) -> bool {
+    #[test]
+    fn inverse_is_identity(i in unit_dual_quaternion(), p in point3(), v in vector3()) {
         let ii = i.inverse();
 
-        relative_eq!(i * ii, UnitDualQuaternion::identity(), epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(i * ii, UnitDualQuaternion::identity(), epsilon = 1.0e-7)
             && relative_eq!(ii * i, UnitDualQuaternion::identity(), epsilon = 1.0e-7)
             && relative_eq!((i * ii) * p, p, epsilon = 1.0e-7)
             && relative_eq!((ii * i) * p, p, epsilon = 1.0e-7)
             && relative_eq!((i * ii) * v, v, epsilon = 1.0e-7)
-            && relative_eq!((ii * i) * v, v, epsilon = 1.0e-7)
+            && relative_eq!((ii * i) * v, v, epsilon = 1.0e-7));
     }
 
     #[cfg_attr(rustfmt, rustfmt_skip)]
+    #[test]
     fn multiply_equals_alga_transform(
-        dq: UnitDualQuaternion<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>
-    ) -> bool {
-        dq * v == dq.transform_vector(&v)
+        dq in unit_dual_quaternion(),
+        v in vector3(),
+        p in point3()
+    ) {
+        prop_assert!(dq * v == dq.transform_vector(&v)
             && dq * p == dq.transform_point(&p)
             && relative_eq!(
                 dq.inverse() * v,
@@ -41,44 +45,46 @@ quickcheck!(
                 dq.inverse() * p,
                 dq.inverse_transform_point(&p),
                 epsilon = 1.0e-7
-            )
+            ));
     }
 
     #[cfg_attr(rustfmt, rustfmt_skip)]
+    #[test]
     fn composition(
-        dq: UnitDualQuaternion<f64>,
-        uq: UnitQuaternion<f64>,
-        t: Translation3<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>
-    ) -> bool {
+        dq in unit_dual_quaternion(),
+        uq in unit_quaternion(),
+        t in translation3(),
+        v in vector3(),
+        p in point3()
+    ) {
         // (rotation × dual quaternion) * point = rotation × (dual quaternion * point)
-        relative_eq!((uq * dq) * v, uq * (dq * v), epsilon = 1.0e-7) &&
-        relative_eq!((uq * dq) * p, uq * (dq * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((uq * dq) * v, uq * (dq * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((uq * dq) * p, uq * (dq * p), epsilon = 1.0e-7));
 
         // (dual quaternion × rotation) * point = dual quaternion × (rotation * point)
-        relative_eq!((dq * uq) * v, dq * (uq * v), epsilon = 1.0e-7) &&
-        relative_eq!((dq * uq) * p, dq * (uq * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((dq * uq) * v, dq * (uq * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((dq * uq) * p, dq * (uq * p), epsilon = 1.0e-7));
 
         // (translation × dual quaternion) * point = translation × (dual quaternion * point)
-        relative_eq!((t * dq) * v,     (dq * v), epsilon = 1.0e-7) &&
-        relative_eq!((t * dq) * p, t * (dq * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((t * dq) * v,     (dq * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((t * dq) * p, t * (dq * p), epsilon = 1.0e-7));
 
         // (dual quaternion × translation) * point = dual quaternion × (translation * point)
-        relative_eq!((dq * t) * v, dq * v,       epsilon = 1.0e-7) &&
-        relative_eq!((dq * t) * p, dq * (t * p), epsilon = 1.0e-7)
+        prop_assert!(relative_eq!((dq * t) * v, dq * v,       epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((dq * t) * p, dq * (t * p), epsilon = 1.0e-7));
     }
 
     #[cfg_attr(rustfmt, rustfmt_skip)]
+    #[test]
     fn all_op_exist(
-        dq: DualQuaternion<f64>,
-        udq: UnitDualQuaternion<f64>,
-        uq: UnitQuaternion<f64>,
-        s: f64,
-        t: Translation3<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>
-    ) -> bool {
+        dq in dual_quaternion(),
+        udq in unit_dual_quaternion(),
+        uq in unit_quaternion(),
+        s in PROPTEST_F64,
+        t in translation3(),
+        v in vector3(),
+        p in point3()
+    ) {
         let dqMs: DualQuaternion<_> = dq * s;
 
         let dqMdq: DualQuaternion<_> = dq * dq;
@@ -145,7 +151,7 @@ quickcheck!(
         iDuq1 /= uq;
         iDuq2 /= &uq;
 
-        dqMs == dqMs1
+        prop_assert!(dqMs == dqMs1
             && dqMdq == dqMdq1
             && dqMdq == dqMdq2
             && dqMudq == dqMudq1
@@ -199,6 +205,6 @@ quickcheck!(
             && uqMi == &uq * udq
             && uqDi == &uq / &udq
             && uqDi == uq / &udq
-            && uqDi == &uq / udq
+            && uqDi == &uq / udq)
     }
 );

tests/geometry/isometry.rs

@@ -1,67 +1,74 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 #![allow(non_snake_case)]
 
-use na::{
-    Isometry2, Isometry3, Point2, Point3, Rotation2, Rotation3, Translation2, Translation3,
-    UnitComplex, UnitQuaternion, Vector2, Vector3,
-};
+use na::{Isometry3, Point3, Vector3};
 
-quickcheck!(
-    fn append_rotation_wrt_point_to_id(r: UnitQuaternion<f64>, p: Point3<f64>) -> bool {
+use crate::proptest::*;
+use proptest::{prop_assert, prop_assert_eq, proptest};
+
+proptest!(
+    #[test]
+    fn append_rotation_wrt_point_to_id(r in unit_quaternion(), p in point3()) {
         let mut iso = Isometry3::identity();
         iso.append_rotation_wrt_point_mut(&r, &p);
 
-        iso == Isometry3::rotation_wrt_point(r, p)
+        prop_assert_eq!(iso, Isometry3::rotation_wrt_point(r, p))
     }
 
-    fn rotation_wrt_point_invariance(r: UnitQuaternion<f64>, p: Point3<f64>) -> bool {
+    #[test]
+    fn rotation_wrt_point_invariance(r in unit_quaternion(), p in point3()) {
         let iso = Isometry3::rotation_wrt_point(r, p);
 
-        relative_eq!(iso * p, p, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(iso * p, p, epsilon = 1.0e-7))
     }
 
-    fn look_at_rh_3(eye: Point3<f64>, target: Point3<f64>, up: Vector3<f64>) -> bool {
+    #[test]
+    fn look_at_rh_3(eye in point3(), target in point3(), up in vector3()) {
         let viewmatrix = Isometry3::look_at_rh(&eye, &target, &up);
-
         let origin = Point3::origin();
-        relative_eq!(viewmatrix * eye, origin, epsilon = 1.0e-7)
+
+        prop_assert!(relative_eq!(viewmatrix * eye, origin, epsilon = 1.0e-7)
             && relative_eq!(
                 (viewmatrix * (target - eye)).normalize(),
                 -Vector3::z(),
                 epsilon = 1.0e-7
-            )
+            ))
     }
 
-    fn observer_frame_3(eye: Point3<f64>, target: Point3<f64>, up: Vector3<f64>) -> bool {
+    #[test]
+    fn observer_frame_3(eye in point3(), target in point3(), up in vector3()) {
         let observer = Isometry3::face_towards(&eye, &target, &up);
-
         let origin = Point3::origin();
-        relative_eq!(observer * origin, eye, epsilon = 1.0e-7)
+
+        prop_assert!(relative_eq!(observer * origin, eye, epsilon = 1.0e-7)
             && relative_eq!(
                 observer * Vector3::z(),
                 (target - eye).normalize(),
                 epsilon = 1.0e-7
-            )
+            ))
     }
 
-    fn inverse_is_identity(i: Isometry3<f64>, p: Point3<f64>, v: Vector3<f64>) -> bool {
+    #[test]
+    fn inverse_is_identity(i in isometry3(), p in point3(), v in vector3()) {
         let ii = i.inverse();
 
-        relative_eq!(i * ii, Isometry3::identity(), epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(i * ii, Isometry3::identity(), epsilon = 1.0e-7)
             && relative_eq!(ii * i, Isometry3::identity(), epsilon = 1.0e-7)
             && relative_eq!((i * ii) * p, p, epsilon = 1.0e-7)
             && relative_eq!((ii * i) * p, p, epsilon = 1.0e-7)
             && relative_eq!((i * ii) * v, v, epsilon = 1.0e-7)
-            && relative_eq!((ii * i) * v, v, epsilon = 1.0e-7)
+            && relative_eq!((ii * i) * v, v, epsilon = 1.0e-7))
     }
 
-    fn inverse_is_parts_inversion(t: Translation3<f64>, r: UnitQuaternion<f64>) -> bool {
+    #[test]
+    fn inverse_is_parts_inversion(t in translation3(), r in unit_quaternion()) {
         let i = t * r;
-        i.inverse() == r.inverse() * t.inverse()
+        prop_assert!(i.inverse() == r.inverse() * t.inverse())
     }
 
-    fn multiply_equals_alga_transform(i: Isometry3<f64>, v: Vector3<f64>, p: Point3<f64>) -> bool {
-        i * v == i.transform_vector(&v)
+    #[test]
+    fn multiply_equals_alga_transform(i in isometry3(), v in vector3(), p in point3()) {
+        prop_assert!(i * v == i.transform_vector(&v)
             && i * p == i.transform_point(&p)
             && relative_eq!(
                 i.inverse() * v,
@@ -72,94 +79,97 @@ quickcheck!(
                 i.inverse() * p,
                 i.inverse_transform_point(&p),
                 epsilon = 1.0e-7
-            )
+            ))
     }
 
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn composition2(
-        i: Isometry2<f64>,
-        uc: UnitComplex<f64>,
-        r: Rotation2<f64>,
-        t: Translation2<f64>,
-        v: Vector2<f64>,
-        p: Point2<f64>
-    ) -> bool {
+        i in isometry2(),
+        uc in unit_complex(),
+        r in rotation2(),
+        t in translation2(),
+        v in vector2(),
+        p in point2()
+    ) {
         // (rotation × translation) * point = rotation × (translation * point)
-        relative_eq!((uc * t) * v, uc * v, epsilon = 1.0e-7)       &&
-        relative_eq!((r  * t) * v, r  * v, epsilon = 1.0e-7)       &&
-        relative_eq!((uc * t) * p, uc * (t * p), epsilon = 1.0e-7) &&
-        relative_eq!((r * t)  * p, r  * (t * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((uc * t) * v, uc * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((r  * t) * v, r  * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((uc * t) * p, uc * (t * p), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((r * t)  * p, r  * (t * p), epsilon = 1.0e-7));
 
         // (translation × rotation) * point = translation × (rotation * point)
-        (t * uc) * v == uc * v       &&
-        (t * r)  * v == r  * v       &&
-        (t * uc) * p == t * (uc * p) &&
-        (t * r)  * p == t * (r  * p) &&
+        prop_assert_eq!((t * uc) * v, uc * v);
+        prop_assert_eq!((t * r)  * v, r  * v);
+        prop_assert_eq!((t * uc) * p, t * (uc * p));
+        prop_assert_eq!((t * r)  * p, t * (r  * p));
 
         // (rotation × isometry) * point = rotation × (isometry * point)
-        relative_eq!((uc * i) * v, uc * (i * v), epsilon = 1.0e-7) &&
-        relative_eq!((uc * i) * p, uc * (i * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((uc * i) * v, uc * (i * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((uc * i) * p, uc * (i * p), epsilon = 1.0e-7));
 
         // (isometry × rotation) * point = isometry × (rotation * point)
-        relative_eq!((i * uc) * v, i * (uc * v), epsilon = 1.0e-7) &&
-        relative_eq!((i * uc) * p, i * (uc * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((i * uc) * v, i * (uc * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((i * uc) * p, i * (uc * p), epsilon = 1.0e-7));
 
         // (translation × isometry) * point = translation × (isometry * point)
-        relative_eq!((t * i) * v,     (i * v), epsilon = 1.0e-7) &&
-        relative_eq!((t * i) * p, t * (i * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((t * i) * v,     (i * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((t * i) * p, t * (i * p), epsilon = 1.0e-7));
 
         // (isometry × translation) * point = isometry × (translation * point)
-        relative_eq!((i * t) * v, i * v,       epsilon = 1.0e-7) &&
-        relative_eq!((i * t) * p, i * (t * p), epsilon = 1.0e-7)
+        prop_assert!(relative_eq!((i * t) * v, i * v,       epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((i * t) * p, i * (t * p), epsilon = 1.0e-7));
     }
 
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn composition3(
-        i: Isometry3<f64>,
-        uq: UnitQuaternion<f64>,
-        r: Rotation3<f64>,
-        t: Translation3<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>
-    ) -> bool {
+        i in isometry3(),
+        uq in unit_quaternion(),
+        r in rotation3(),
+        t in translation3(),
+        v in vector3(),
+        p in point3()
+    ) {
         // (rotation × translation) * point = rotation × (translation * point)
-        relative_eq!((uq * t) * v, uq * v, epsilon = 1.0e-7)       &&
-        relative_eq!((r  * t) * v, r  * v, epsilon = 1.0e-7)       &&
-        relative_eq!((uq * t) * p, uq * (t * p), epsilon = 1.0e-7) &&
-        relative_eq!((r * t)  * p, r  * (t * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((uq * t) * v, uq * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((r  * t) * v, r  * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((uq * t) * p, uq * (t * p), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((r * t)  * p, r  * (t * p), epsilon = 1.0e-7));
 
         // (translation × rotation) * point = translation × (rotation * point)
-        (t * uq) * v == uq * v       &&
-        (t * r)  * v == r  * v       &&
-        (t * uq) * p == t * (uq * p) &&
-        (t * r)  * p == t * (r  * p) &&
+        prop_assert_eq!((t * uq) * v, uq * v);
+        prop_assert_eq!((t * r)  * v, r  * v);
+        prop_assert_eq!((t * uq) * p, t * (uq * p));
+        prop_assert_eq!((t * r)  * p, t * (r  * p));
 
         // (rotation × isometry) * point = rotation × (isometry * point)
-        relative_eq!((uq * i) * v, uq * (i * v), epsilon = 1.0e-7) &&
-        relative_eq!((uq * i) * p, uq * (i * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((uq * i) * v, uq * (i * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((uq * i) * p, uq * (i * p), epsilon = 1.0e-7));
 
         // (isometry × rotation) * point = isometry × (rotation * point)
-        relative_eq!((i * uq) * v, i * (uq * v), epsilon = 1.0e-7) &&
-        relative_eq!((i * uq) * p, i * (uq * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((i * uq) * v, i * (uq * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((i * uq) * p, i * (uq * p), epsilon = 1.0e-7));
 
         // (translation × isometry) * point = translation × (isometry * point)
-        relative_eq!((t * i) * v,     (i * v), epsilon = 1.0e-7) &&
-        relative_eq!((t * i) * p, t * (i * p), epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!((t * i) * v,     (i * v), epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((t * i) * p, t * (i * p), epsilon = 1.0e-7));
 
         // (isometry × translation) * point = isometry × (translation * point)
-        relative_eq!((i * t) * v, i * v,       epsilon = 1.0e-7) &&
-        relative_eq!((i * t) * p, i * (t * p), epsilon = 1.0e-7)
+        prop_assert!(relative_eq!((i * t) * v, i * v,       epsilon = 1.0e-7));
+        prop_assert!(relative_eq!((i * t) * p, i * (t * p), epsilon = 1.0e-7));
     }
 
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn all_op_exist(
-        i: Isometry3<f64>,
-        uq: UnitQuaternion<f64>,
-        t: Translation3<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>,
-        r: Rotation3<f64>
-    ) -> bool {
+        i in isometry3(),
+        uq in unit_quaternion(),
+        t in translation3(),
+        v in vector3(),
+        p in point3(),
+        r in rotation3()
+    ) {
         let iMi = i * i;
         let iMuq = i * uq;
         let iDi = i / i;
@@ -210,7 +220,7 @@ quickcheck!(
         iDuq1 /= uq;
         iDuq2 /= &uq;
 
-        iMt == iMt1
+        prop_assert!(iMt == iMt1
             && iMt == iMt2
             && iMi == iMi1
             && iMi == iMi2
@@ -261,6 +271,6 @@ quickcheck!(
             && rMt == &r * t
             && uqMt == &uq * &t
             && uqMt == uq * &t
-            && uqMt == &uq * t
+            && uqMt == &uq * t)
     }
 );

tests/geometry/point.rs

@@ -92,11 +92,3 @@ fn to_homogeneous() {
 
     assert_eq!(a.to_homogeneous(), expected);
 }
-
-#[cfg(feature = "arbitrary")]
-quickcheck!(
-    fn point_sub(pt1: Point3<f64>, pt2: Point3<f64>) -> bool {
-        let dpt = &pt2 - &pt1;
-        relative_eq!(pt2, pt1 + dpt, epsilon = 1.0e-7)
-    }
-);

tests/geometry/projection.rs

@@ -33,27 +33,32 @@ fn perspective_matrix_point_transformation() {
     );
 }
 
-#[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
-    use na::{Orthographic3, Perspective3, Point3};
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
+    use na::{Orthographic3, Perspective3};
 
-    quickcheck! {
-        fn perspective_project_unproject(pt: Point3<f64>) -> bool {
+    use crate::proptest::*;
+    use proptest::{prop_assert, proptest};
+
+    proptest! {
+        #[test]
+        fn perspective_project_unproject(pt in point3()) {
             let proj = Perspective3::new(800.0 / 600.0, 3.14 / 2.0, 1.0, 1000.0);
 
             let projected   = proj.project_point(&pt);
             let unprojected = proj.unproject_point(&projected);
 
-            relative_eq!(pt, unprojected, epsilon = 1.0e-7)
+            prop_assert!(relative_eq!(pt, unprojected, epsilon = 1.0e-7))
         }
 
-        fn orthographic_project_unproject(pt: Point3<f64>) -> bool {
+        #[test]
+        fn orthographic_project_unproject(pt in point3()) {
             let proj = Orthographic3::new(1.0, 2.0, -3.0, -2.5, 10.0, 900.0);
 
             let projected   = proj.project_point(&pt);
             let unprojected = proj.unproject_point(&projected);
 
-            relative_eq!(pt, unprojected, epsilon = 1.0e-7)
+            prop_assert!(relative_eq!(pt, unprojected, epsilon = 1.0e-7))
         }
     }
 }

tests/geometry/quaternion.rs

@@ -1,15 +1,19 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 #![allow(non_snake_case)]
 
-use na::{Point3, Quaternion, Rotation3, Unit, UnitQuaternion, Vector3};
+use na::{Unit, UnitQuaternion};
 
-quickcheck!(
+use crate::proptest::*;
+use proptest::{prop_assert, proptest};
+
+proptest!(
     /*
      *
      * Euler angles.
      *
      */
-    fn from_euler_angles(r: f64, p: f64, y: f64) -> bool {
+    #[test]
+    fn from_euler_angles(r in PROPTEST_F64, p in PROPTEST_F64, y in PROPTEST_F64) {
         let roll = UnitQuaternion::from_euler_angles(r, 0.0, 0.0);
         let pitch = UnitQuaternion::from_euler_angles(0.0, p, 0.0);
         let yaw = UnitQuaternion::from_euler_angles(0.0, 0.0, y);
@@ -20,20 +24,21 @@ quickcheck!(
         let rpitch = pitch.to_rotation_matrix();
         let ryaw = yaw.to_rotation_matrix();
 
-        relative_eq!(rroll[(0, 0)],  1.0, epsilon = 1.0e-7) && // rotation wrt. x axis.
-        relative_eq!(rpitch[(1, 1)], 1.0, epsilon = 1.0e-7) && // rotation wrt. y axis.
-        relative_eq!(ryaw[(2, 2)],   1.0, epsilon = 1.0e-7) && // rotation wrt. z axis.
-        relative_eq!(yaw * pitch * roll, rpy, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(rroll[(0, 0)],  1.0, epsilon = 1.0e-7)); // rotation wrt. x axis.
+        prop_assert!(relative_eq!(rpitch[(1, 1)], 1.0, epsilon = 1.0e-7)); // rotation wrt. y axis.
+        prop_assert!(relative_eq!(ryaw[(2, 2)],   1.0, epsilon = 1.0e-7)); // rotation wrt. z axis.
+        prop_assert!(relative_eq!(yaw * pitch * roll, rpy, epsilon = 1.0e-7));
     }
 
-    fn euler_angles(r: f64, p: f64, y: f64) -> bool {
+    #[test]
+    fn euler_angles(r in PROPTEST_F64, p in PROPTEST_F64, y in PROPTEST_F64) {
         let rpy = UnitQuaternion::from_euler_angles(r, p, y);
         let (roll, pitch, yaw) = rpy.euler_angles();
-        relative_eq!(
+        prop_assert!(relative_eq!(
             UnitQuaternion::from_euler_angles(roll, pitch, yaw),
             rpy,
             epsilon = 1.0e-7
-        )
+        ))
     }
 
     /*
@@ -41,12 +46,13 @@ quickcheck!(
      * From/to rotation matrix.
      *
      */
-    fn unit_quaternion_rotation_conversion(q: UnitQuaternion<f64>) -> bool {
+    #[test]
+    fn unit_quaternion_rotation_conversion(q in unit_quaternion()) {
         let r = q.to_rotation_matrix();
         let qq = UnitQuaternion::from_rotation_matrix(&r);
         let rr = qq.to_rotation_matrix();
 
-        relative_eq!(q, qq, epsilon = 1.0e-7) && relative_eq!(r, rr, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(q, qq, epsilon = 1.0e-7) && relative_eq!(r, rr, epsilon = 1.0e-7))
     }
 
     /*
@@ -55,24 +61,25 @@ quickcheck!(
      *
      */
 
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn unit_quaternion_transformation(
-        q: UnitQuaternion<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>
-    ) -> bool {
+        q in unit_quaternion(),
+        v in vector3(),
+        p in point3()
+    ) {
         let r = q.to_rotation_matrix();
         let rv = r * v;
         let rp = r * p;
 
-        relative_eq!(q * v, rv, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(q * v, rv, epsilon = 1.0e-7)
             && relative_eq!(q * &v, rv, epsilon = 1.0e-7)
             && relative_eq!(&q * v, rv, epsilon = 1.0e-7)
             && relative_eq!(&q * &v, rv, epsilon = 1.0e-7)
             && relative_eq!(q * p, rp, epsilon = 1.0e-7)
             && relative_eq!(q * &p, rp, epsilon = 1.0e-7)
             && relative_eq!(&q * p, rp, epsilon = 1.0e-7)
-            && relative_eq!(&q * &p, rp, epsilon = 1.0e-7)
+            && relative_eq!(&q * &p, rp, epsilon = 1.0e-7))
     }
 
     /*
@@ -80,16 +87,17 @@ quickcheck!(
      * Inversion.
      *
      */
-    fn unit_quaternion_inv(q: UnitQuaternion<f64>) -> bool {
+    #[test]
+    fn unit_quaternion_inv(q in unit_quaternion()) {
         let iq = q.inverse();
-        relative_eq!(&iq * &q, UnitQuaternion::identity(), epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(&iq * &q, UnitQuaternion::identity(), epsilon = 1.0e-7)
             && relative_eq!(iq * &q, UnitQuaternion::identity(), epsilon = 1.0e-7)
             && relative_eq!(&iq * q, UnitQuaternion::identity(), epsilon = 1.0e-7)
             && relative_eq!(iq * q, UnitQuaternion::identity(), epsilon = 1.0e-7)
             && relative_eq!(&q * &iq, UnitQuaternion::identity(), epsilon = 1.0e-7)
             && relative_eq!(q * &iq, UnitQuaternion::identity(), epsilon = 1.0e-7)
             && relative_eq!(&q * iq, UnitQuaternion::identity(), epsilon = 1.0e-7)
-            && relative_eq!(q * iq, UnitQuaternion::identity(), epsilon = 1.0e-7)
+            && relative_eq!(q * iq, UnitQuaternion::identity(), epsilon = 1.0e-7))
     }
 
     /*
@@ -97,14 +105,15 @@ quickcheck!(
      * Quaterion * Vector == Rotation * Vector
      *
      */
-    fn unit_quaternion_mul_vector(q: UnitQuaternion<f64>, v: Vector3<f64>, p: Point3<f64>) -> bool {
+    #[test]
+    fn unit_quaternion_mul_vector(q in unit_quaternion(), v in vector3(), p in point3()) {
         let r = q.to_rotation_matrix();
 
-        relative_eq!(q * v, r * v, epsilon = 1.0e-7) &&
-        relative_eq!(q * p, r * p, epsilon = 1.0e-7) &&
+        prop_assert!(relative_eq!(q * v, r * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(q * p, r * p, epsilon = 1.0e-7));
         // Equivalence q = -q
-        relative_eq!(UnitQuaternion::new_unchecked(-q.into_inner()) * v, r * v, epsilon = 1.0e-7) &&
-        relative_eq!(UnitQuaternion::new_unchecked(-q.into_inner()) * p, r * p, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(UnitQuaternion::new_unchecked(-q.into_inner()) * v, r * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(UnitQuaternion::new_unchecked(-q.into_inner()) * p, r * p, epsilon = 1.0e-7));
     }
 
     /*
@@ -112,23 +121,25 @@ quickcheck!(
      * Unit quaternion double-covering.
      *
      */
-    fn unit_quaternion_double_covering(q: UnitQuaternion<f64>) -> bool {
+    #[test]
+    fn unit_quaternion_double_covering(q in unit_quaternion()) {
         let mq = UnitQuaternion::new_unchecked(-q.into_inner());
-        mq == q && mq.angle() == q.angle() && mq.axis() == q.axis()
+        prop_assert!(mq == q && mq.angle() == q.angle() && mq.axis() == q.axis())
     }
 
     // Test that all operators (incl. all combinations of references) work.
     // See the top comment on `geometry/quaternion_ops.rs` for details on which operations are
     // supported.
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn all_op_exist(
-        q: Quaternion<f64>,
-        uq: UnitQuaternion<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>,
-        r: Rotation3<f64>,
-        s: f64
-    ) -> bool {
+        q in quaternion(),
+        uq in unit_quaternion(),
+        v in vector3(),
+        p in point3(),
+        r in rotation3(),
+        s in PROPTEST_F64
+    ) {
         let uv = Unit::new_normalize(v);
 
         let qpq = q + q;
@@ -196,7 +207,7 @@ quickcheck!(
         uqDr1 /= r;
         uqDr2 /= &r;
 
-        qMs1 == qMs
+        prop_assert!(qMs1 == qMs
             && qMq1 == qMq
             && qMq1 == qMq2
             && qpq1 == qpq
@@ -250,6 +261,6 @@ quickcheck!(
             && uqMv == &uq * v
             && uqMuv == &uq * &uv
             && uqMuv == uq * &uv
-            && uqMuv == &uq * uv
+            && uqMuv == &uq * uv)
     }
 );

tests/geometry/rotation.rs

@@ -30,44 +30,50 @@ fn quaternion_euler_angles_issue_494() {
     assert_eq!(angs.2, 0.0);
 }
 
-#[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
-    use na::{self, Rotation2, Rotation3, Unit, Vector2, Vector3};
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
+    use na::{self, Rotation2, Rotation3, Unit};
     use simba::scalar::RealField;
     use std::f64;
 
-    quickcheck! {
+    use crate::proptest::*;
+    use proptest::{prop_assert, prop_assert_eq, proptest};
+
+    proptest! {
         /*
          *
          * Euler angles.
          *
          */
-        fn from_euler_angles(r: f64, p: f64, y: f64) -> bool {
+        #[test]
+        fn from_euler_angles(r in PROPTEST_F64, p in PROPTEST_F64, y in PROPTEST_F64) {
             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);
 
             let rpy = Rotation3::from_euler_angles(r, p, 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
+            prop_assert_eq!(roll[(0, 0)], 1.0); // rotation wrt. x axis.
+            prop_assert_eq!(pitch[(1, 1)], 1.0); // rotation wrt. y axis.
+            prop_assert_eq!(yaw[(2, 2)], 1.0); // rotation wrt. z axis.
+            prop_assert_eq!(yaw * pitch * roll, rpy);
         }
 
-        fn euler_angles(r: f64, p: f64, y: f64) -> bool {
+        #[test]
+        fn euler_angles(r in PROPTEST_F64, p in PROPTEST_F64, y in PROPTEST_F64) {
             let rpy = Rotation3::from_euler_angles(r, p, y);
             let (roll, pitch, yaw) = rpy.euler_angles();
-            relative_eq!(Rotation3::from_euler_angles(roll, pitch, yaw), rpy, epsilon = 1.0e-7)
+            prop_assert!(relative_eq!(Rotation3::from_euler_angles(roll, pitch, yaw), rpy, epsilon = 1.0e-7));
         }
 
-        fn euler_angles_gimble_lock(r: f64, y: f64) -> bool {
+        #[test]
+        fn euler_angles_gimble_lock(r in PROPTEST_F64, y in PROPTEST_F64) {
             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.euler_angles();
             let (neg_r, neg_p, neg_y) = neg.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)
+            prop_assert!(relative_eq!(Rotation3::from_euler_angles(pos_r, pos_p, pos_y), pos, epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(Rotation3::from_euler_angles(neg_r, neg_p, neg_y), neg, epsilon = 1.0e-7));
         }
 
         /*
@@ -75,26 +81,28 @@ mod quickcheck_tests {
          * Inversion is transposition.
          *
          */
-        fn rotation_inv_3(a: Rotation3<f64>) -> bool {
+        #[test]
+        fn rotation_inv_3(a in rotation3()) {
             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)
+            prop_assert_eq!(ta, ia);
+            prop_assert!(relative_eq!(&ta * &a,  Rotation3::identity(), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(&ia *  a,  Rotation3::identity(), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(  a * &ta, Rotation3::identity(), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(  a *  ia, Rotation3::identity(), epsilon = 1.0e-7));
         }
 
-        fn rotation_inv_2(a: Rotation2<f64>) -> bool {
+        #[test]
+        fn rotation_inv_2(a in rotation2()) {
             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)
+            prop_assert_eq!(ta, ia);
+            prop_assert!(relative_eq!(&ta * &a,  Rotation2::identity(), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(&ia *  a,  Rotation2::identity(), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(  a * &ta, Rotation2::identity(), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(  a *  ia, Rotation2::identity(), epsilon = 1.0e-7));
         }
 
         /*
@@ -102,12 +110,14 @@ mod quickcheck_tests {
          * Angle between vectors.
          *
          */
-        fn angle_is_commutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
-            a.angle(&b) == b.angle(&a)
+        #[test]
+        fn angle_is_commutative_2(a in vector2(), b in vector2()) {
+            prop_assert_eq!(a.angle(&b), b.angle(&a))
         }
 
-        fn angle_is_commutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
-            a.angle(&b) == b.angle(&a)
+        #[test]
+        fn angle_is_commutative_3(a in vector3(), b in vector3()) {
+            prop_assert_eq!(a.angle(&b), b.angle(&a))
         }
 
         /*
@@ -115,50 +125,46 @@ mod quickcheck_tests {
          * Rotation matrix between vectors.
          *
          */
-        fn rotation_between_is_anticommutative_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
+        #[test]
+        fn rotation_between_is_anticommutative_2(a in vector2(), b in vector2()) {
             let rab = Rotation2::rotation_between(&a, &b);
             let rba = Rotation2::rotation_between(&b, &a);
 
-            relative_eq!(rab * rba, Rotation2::identity())
+            prop_assert!(relative_eq!(rab * rba, Rotation2::identity()));
         }
 
-        fn rotation_between_is_anticommutative_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
+        #[test]
+        fn rotation_between_is_anticommutative_3(a in vector3(), b in vector3()) {
             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
+                prop_assert!(relative_eq!(rab * rba, Rotation3::identity(), epsilon = 1.0e-7));
             }
         }
 
-        fn rotation_between_is_identity(v2: Vector2<f64>, v3: Vector3<f64>) -> bool {
+        #[test]
+        fn rotation_between_is_identity(v2 in vector2(), v3 in vector3()) {
             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()
+            prop_assert!(relative_eq!(v2.angle(&vv2), 0.0, epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(v3.angle(&vv3), 0.0, epsilon = 1.0e-7));
+            prop_assert!(relative_eq!(Rotation2::rotation_between(&v2, &vv2), Rotation2::identity()));
+            prop_assert_eq!(Rotation3::rotation_between(&v3, &vv3).unwrap(), Rotation3::identity());
         }
 
-        fn rotation_between_2(a: Vector2<f64>, b: Vector2<f64>) -> bool {
+        #[test]
+        fn rotation_between_2(a in vector2(), b in vector2()) {
             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
+                prop_assert!(relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7))
             }
         }
 
-        fn rotation_between_3(a: Vector3<f64>, b: Vector3<f64>) -> bool {
+        #[test]
+        fn rotation_between_3(a in vector3(), b in vector3()) {
             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
+                prop_assert!(relative_eq!((r * a).angle(&b), 0.0, epsilon = 1.0e-7))
             }
         }
 
@@ -168,25 +174,27 @@ mod quickcheck_tests {
          * Rotation construction.
          *
          */
-        fn new_rotation_2(angle: f64) -> bool {
+        #[test]
+        fn new_rotation_2(angle in PROPTEST_F64) {
             let r = Rotation2::new(angle);
 
             let angle = na::wrap(angle, -f64::pi(), f64::pi());
-            relative_eq!(r.angle(), angle, epsilon = 1.0e-7)
+            prop_assert!(relative_eq!(r.angle(), angle, epsilon = 1.0e-7))
         }
 
-        fn new_rotation_3(axisangle: Vector3<f64>) -> bool {
+        #[test]
+        fn new_rotation_3(axisangle in vector3()) {
             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) &&
+                prop_assert!((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.axis().unwrap(), -axis, epsilon = 1.0e-7)))
             }
             else {
-                r == Rotation3::identity()
+                prop_assert_eq!(r, Rotation3::identity())
             }
         }
 
@@ -195,28 +203,30 @@ mod quickcheck_tests {
          * Rotation pow.
          *
          */
-        fn powf_rotation_2(angle: f64, pow: f64) -> bool {
+        #[test]
+        fn powf_rotation_2(angle in PROPTEST_F64, pow in PROPTEST_F64) {
             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)
+            prop_assert!(relative_eq!(r.angle(), pangle, epsilon = 1.0e-7));
         }
 
-        fn powf_rotation_3(axisangle: Vector3<f64>, pow: f64) -> bool {
+        #[test]
+        fn powf_rotation_3(axisangle in vector3(), pow in PROPTEST_F64) {
             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) &&
+                prop_assert!((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.axis().unwrap(), -axis, epsilon = 1.0e-7)));
             }
             else {
-                r == Rotation3::identity()
+                prop_assert_eq!(r, Rotation3::identity())
             }
         }
     }

tests/geometry/similarity.rs

@@ -1,41 +1,45 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 #![allow(non_snake_case)]
 
-use na::{Isometry3, Point3, Similarity3, Translation3, UnitQuaternion, Vector3};
+use na::Similarity3;
 
-quickcheck!(
-    fn inverse_is_identity(i: Similarity3<f64>, p: Point3<f64>, v: Vector3<f64>) -> bool {
+use crate::proptest::*;
+use proptest::{prop_assert, prop_assert_eq, proptest};
+
+proptest!(
+    #[test]
+    fn inverse_is_identity(i in similarity3(), p in point3(), v in vector3()) {
         let ii = i.inverse();
 
-        relative_eq!(i * ii, Similarity3::identity(), epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(i * ii, Similarity3::identity(), epsilon = 1.0e-7)
             && relative_eq!(ii * i, Similarity3::identity(), epsilon = 1.0e-7)
             && relative_eq!((i * ii) * p, p, epsilon = 1.0e-7)
             && relative_eq!((ii * i) * p, p, epsilon = 1.0e-7)
             && relative_eq!((i * ii) * v, v, epsilon = 1.0e-7)
-            && relative_eq!((ii * i) * v, v, epsilon = 1.0e-7)
+            && relative_eq!((ii * i) * v, v, epsilon = 1.0e-7))
     }
 
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn inverse_is_parts_inversion(
-        t: Translation3<f64>,
-        r: UnitQuaternion<f64>,
-        scaling: f64
-    ) -> bool {
-        if relative_eq!(scaling, 0.0) {
-            true
-        } else {
+        t in translation3(),
+        r in unit_quaternion(),
+        scaling in PROPTEST_F64
+    ) {
+        if !relative_eq!(scaling, 0.0) {
             let s = Similarity3::from_isometry(t * r, scaling);
-            s.inverse() == Similarity3::from_scaling(1.0 / scaling) * r.inverse() * t.inverse()
+            prop_assert_eq!(s.inverse(), Similarity3::from_scaling(1.0 / scaling) * r.inverse() * t.inverse())
         }
     }
 
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn multiply_equals_alga_transform(
-        s: Similarity3<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>
-    ) -> bool {
-        s * v == s.transform_vector(&v)
+        s in similarity3(),
+        v in vector3(),
+        p in point3()
+    ) {
+        prop_assert!(s * v == s.transform_vector(&v)
             && s * p == s.transform_point(&p)
             && relative_eq!(
                 s.inverse() * v,
@@ -46,114 +50,114 @@ quickcheck!(
                 s.inverse() * p,
                 s.inverse_transform_point(&p),
                 epsilon = 1.0e-7
-            )
+            ))
     }
 
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn composition(
-        i: Isometry3<f64>,
-        uq: UnitQuaternion<f64>,
-        t: Translation3<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>,
-        scaling: f64
-    ) -> bool {
-        if relative_eq!(scaling, 0.0) {
-            return true;
-        }
-
+        i in isometry3(),
+        uq in unit_quaternion(),
+        t in translation3(),
+        v in vector3(),
+        p in point3(),
+        scaling in PROPTEST_F64
+    ) {
+        if !relative_eq!(scaling, 0.0) {
             let s = Similarity3::from_scaling(scaling);
 
             // (rotation × translation × scaling) × point = rotation × (translation × (scaling × point))
-        relative_eq!((uq * t * s) * v, uq * (scaling * v), epsilon = 1.0e-7)       &&
-        relative_eq!((uq * t * s) * p, uq * (t * (scaling * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((uq * t * s) * v, uq * (scaling * v), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((uq * t * s) * p, uq * (t * (scaling * p)), epsilon = 1.0e-7));
 
             // (translation × rotation × scaling) × point = translation × (rotation × (scaling × point))
-        relative_eq!((t * uq * s) * v, uq * (scaling * v), epsilon = 1.0e-7)       &&
-        relative_eq!((t * uq * s) * p, t * (uq * (scaling * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((t * uq * s) * v, uq * (scaling * v), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((t * uq * s) * p, t * (uq * (scaling * p)), epsilon = 1.0e-7));
 
             // (rotation × isometry × scaling) × point = rotation × (isometry × (scaling × point))
-        relative_eq!((uq * i * s) * v, uq * (i * (scaling * v)), epsilon = 1.0e-7) &&
-        relative_eq!((uq * i * s) * p, uq * (i * (scaling * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((uq * i * s) * v, uq * (i * (scaling * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((uq * i * s) * p, uq * (i * (scaling * p)), epsilon = 1.0e-7));
 
             // (isometry × rotation × scaling) × point = isometry × (rotation × (scaling × point))
-        relative_eq!((i * uq * s) * v, i * (uq * (scaling * v)), epsilon = 1.0e-7) &&
-        relative_eq!((i * uq * s) * p, i * (uq * (scaling * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((i * uq * s) * v, i * (uq * (scaling * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((i * uq * s) * p, i * (uq * (scaling * p)), epsilon = 1.0e-7));
 
             // (translation × isometry × scaling) × point = translation × (isometry × (scaling × point))
-        relative_eq!((t * i * s) * v,     (i * (scaling * v)), epsilon = 1.0e-7) &&
-        relative_eq!((t * i * s) * p, t * (i * (scaling * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((t * i * s) * v,     (i * (scaling * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((t * i * s) * p, t * (i * (scaling * p)), epsilon = 1.0e-7));
 
             // (isometry × translation × scaling) × point = isometry × (translation × (scaling × point))
-        relative_eq!((i * t * s) * v, i * (scaling * v),       epsilon = 1.0e-7) &&
-        relative_eq!((i * t * s) * p, i * (t * (scaling * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((i * t * s) * v, i * (scaling * v),       epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((i * t * s) * p, i * (t * (scaling * p)), epsilon = 1.0e-7));
 
 
             /*
              * Same as before but with scaling on the middle.
              */
             // (rotation × scaling × translation) × point = rotation × (scaling × (translation × point))
-        relative_eq!((uq * s * t) * v, uq * (scaling * v), epsilon = 1.0e-7)       &&
-        relative_eq!((uq * s * t) * p, uq * (scaling * (t * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((uq * s * t) * v, uq * (scaling * v), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((uq * s * t) * p, uq * (scaling * (t * p)), epsilon = 1.0e-7));
 
             // (translation × scaling × rotation) × point = translation × (scaling × (rotation × point))
-        relative_eq!((t * s * uq) * v, scaling * (uq * v), epsilon = 1.0e-7)       &&
-        relative_eq!((t * s * uq) * p, t * (scaling * (uq * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((t * s * uq) * v, scaling * (uq * v), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((t * s * uq) * p, t * (scaling * (uq * p)), epsilon = 1.0e-7));
 
             // (rotation × scaling × isometry) × point = rotation × (scaling × (isometry × point))
-        relative_eq!((uq * s * i) * v, uq * (scaling * (i * v)), epsilon = 1.0e-7) &&
-        relative_eq!((uq * s * i) * p, uq * (scaling * (i * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((uq * s * i) * v, uq * (scaling * (i * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((uq * s * i) * p, uq * (scaling * (i * p)), epsilon = 1.0e-7));
 
             // (isometry × scaling × rotation) × point = isometry × (scaling × (rotation × point))
-        relative_eq!((i * s * uq) * v, i * (scaling * (uq * v)), epsilon = 1.0e-7) &&
-        relative_eq!((i * s * uq) * p, i * (scaling * (uq * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((i * s * uq) * v, i * (scaling * (uq * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((i * s * uq) * p, i * (scaling * (uq * p)), epsilon = 1.0e-7));
 
             // (translation × scaling × isometry) × point = translation × (scaling × (isometry × point))
-        relative_eq!((t * s * i) * v,     (scaling * (i * v)), epsilon = 1.0e-7) &&
-        relative_eq!((t * s * i) * p, t * (scaling * (i * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((t * s * i) * v,     (scaling * (i * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((t * s * i) * p, t * (scaling * (i * p)), epsilon = 1.0e-7));
 
             // (isometry × scaling × translation) × point = isometry × (scaling × (translation × point))
-        relative_eq!((i * s * t) * v, i * (scaling * v),       epsilon = 1.0e-7) &&
-        relative_eq!((i * s * t) * p, i * (scaling * (t * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((i * s * t) * v, i * (scaling * v),       epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((i * s * t) * p, i * (scaling * (t * p)), epsilon = 1.0e-7));
 
 
             /*
              * Same as before but with scaling on the left.
              */
             // (scaling × rotation × translation) × point = scaling × (rotation × (translation × point))
-        relative_eq!((s * uq * t) * v, scaling * (uq * v), epsilon = 1.0e-7)       &&
-        relative_eq!((s * uq * t) * p, scaling * (uq * (t * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((s * uq * t) * v, scaling * (uq * v), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((s * uq * t) * p, scaling * (uq * (t * p)), epsilon = 1.0e-7));
 
             // (scaling × translation × rotation) × point = scaling × (translation × (rotation × point))
-        relative_eq!((s * t * uq) * v, scaling * (uq * v), epsilon = 1.0e-7)       &&
-        relative_eq!((s * t * uq) * p, scaling * (t * (uq * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((s * t * uq) * v, scaling * (uq * v), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((s * t * uq) * p, scaling * (t * (uq * p)), epsilon = 1.0e-7));
 
             // (scaling × rotation × isometry) × point = scaling × (rotation × (isometry × point))
-        relative_eq!((s * uq * i) * v, scaling * (uq * (i * v)), epsilon = 1.0e-7) &&
-        relative_eq!((s * uq * i) * p, scaling * (uq * (i * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((s * uq * i) * v, scaling * (uq * (i * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((s * uq * i) * p, scaling * (uq * (i * p)), epsilon = 1.0e-7));
 
             // (scaling × isometry × rotation) × point = scaling × (isometry × (rotation × point))
-        relative_eq!((s * i * uq) * v, scaling * (i * (uq * v)), epsilon = 1.0e-7) &&
-        relative_eq!((s * i * uq) * p, scaling * (i * (uq * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((s * i * uq) * v, scaling * (i * (uq * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((s * i * uq) * p, scaling * (i * (uq * p)), epsilon = 1.0e-7));
 
             // (scaling × translation × isometry) × point = scaling × (translation × (isometry × point))
-        relative_eq!((s * t * i) * v,     (scaling * (i * v)), epsilon = 1.0e-7) &&
-        relative_eq!((s * t * i) * p, scaling * (t * (i * p)), epsilon = 1.0e-7) &&
+            prop_assert!(relative_eq!((s * t * i) * v,     (scaling * (i * v)), epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((s * t * i) * p, scaling * (t * (i * p)), epsilon = 1.0e-7));
 
             // (scaling × isometry × translation) × point = scaling × (isometry × (translation × point))
-        relative_eq!((s * i * t) * v, scaling * (i * v),       epsilon = 1.0e-7) &&
-        relative_eq!((s * i * t) * p, scaling * (i * (t * p)), epsilon = 1.0e-7)
+            prop_assert!(relative_eq!((s * i * t) * v, scaling * (i * v),       epsilon = 1.0e-7));
+            prop_assert!(relative_eq!((s * i * t) * p, scaling * (i * (t * p)), epsilon = 1.0e-7));
+        }
     }
 
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn all_op_exist(
-        s: Similarity3<f64>,
-        i: Isometry3<f64>,
-        uq: UnitQuaternion<f64>,
-        t: Translation3<f64>,
-        v: Vector3<f64>,
-        p: Point3<f64>
-    ) -> bool {
+        s in similarity3(),
+        i in isometry3(),
+        uq in unit_quaternion(),
+        t in translation3(),
+        v in vector3(),
+        p in point3()
+    ) {
         let sMs = s * s;
         let sMuq = s * uq;
         let sDs = s / s;
@@ -216,7 +220,7 @@ quickcheck!(
         sDi1 /= i;
         sDi2 /= &i;
 
-        sMt == sMt1
+        prop_assert!(sMt == sMt1
             && sMt == sMt2
             && sMs == sMs1
             && sMs == sMs2
@@ -271,6 +275,6 @@ quickcheck!(
             && iMs == &i * s
             && iDs == &i / &s
             && iDs == i / &s
-            && iDs == &i / s
+            && iDs == &i / s)
     }
 );

tests/geometry/unit_complex.rs

@@ -1,20 +1,25 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 #![allow(non_snake_case)]
 
-use na::{Point2, Rotation2, Unit, UnitComplex, Vector2};
+use na::{Unit, UnitComplex};
 
-quickcheck!(
+use crate::proptest::*;
+use proptest::{prop_assert, proptest};
+
+proptest!(
     /*
      *
      * From/to rotation matrix.
      *
      */
-    fn unit_complex_rotation_conversion(c: UnitComplex<f64>) -> bool {
+    #[test]
+    fn unit_complex_rotation_conversion(c in unit_complex()) {
         let r = c.to_rotation_matrix();
         let cc = UnitComplex::from_rotation_matrix(&r);
         let rr = cc.to_rotation_matrix();
 
-        relative_eq!(c, cc, epsilon = 1.0e-7) && relative_eq!(r, rr, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(c, cc, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(r, rr, epsilon = 1.0e-7));
     }
 
     /*
@@ -22,19 +27,20 @@ quickcheck!(
      * Point/Vector transformation.
      *
      */
-    fn unit_complex_transformation(c: UnitComplex<f64>, v: Vector2<f64>, p: Point2<f64>) -> bool {
+    #[test]
+    fn unit_complex_transformation(c in unit_complex(), v in vector2(), p in point2()) {
         let r = c.to_rotation_matrix();
         let rv = r * v;
         let rp = r * p;
 
-        relative_eq!(c * v, rv, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(c * v, rv, epsilon = 1.0e-7)
             && relative_eq!(c * &v, rv, epsilon = 1.0e-7)
             && relative_eq!(&c * v, rv, epsilon = 1.0e-7)
             && relative_eq!(&c * &v, rv, epsilon = 1.0e-7)
             && relative_eq!(c * p, rp, epsilon = 1.0e-7)
             && relative_eq!(c * &p, rp, epsilon = 1.0e-7)
             && relative_eq!(&c * p, rp, epsilon = 1.0e-7)
-            && relative_eq!(&c * &p, rp, epsilon = 1.0e-7)
+            && relative_eq!(&c * &p, rp, epsilon = 1.0e-7))
     }
 
     /*
@@ -42,39 +48,43 @@ quickcheck!(
      * Inversion.
      *
      */
-    fn unit_complex_inv(c: UnitComplex<f64>) -> bool {
+    #[test]
+    fn unit_complex_inv(c in unit_complex()) {
         let iq = c.inverse();
-        relative_eq!(&iq * &c, UnitComplex::identity(), epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(&iq * &c, UnitComplex::identity(), epsilon = 1.0e-7)
             && relative_eq!(iq * &c, UnitComplex::identity(), epsilon = 1.0e-7)
             && relative_eq!(&iq * c, UnitComplex::identity(), epsilon = 1.0e-7)
             && relative_eq!(iq * c, UnitComplex::identity(), epsilon = 1.0e-7)
             && relative_eq!(&c * &iq, UnitComplex::identity(), epsilon = 1.0e-7)
             && relative_eq!(c * &iq, UnitComplex::identity(), epsilon = 1.0e-7)
             && relative_eq!(&c * iq, UnitComplex::identity(), epsilon = 1.0e-7)
-            && relative_eq!(c * iq, UnitComplex::identity(), epsilon = 1.0e-7)
+            && relative_eq!(c * iq, UnitComplex::identity(), epsilon = 1.0e-7))
     }
 
     /*
      *
-     * Quaterion * Vector == Rotation * Vector
+     * Quaternion * Vector == Rotation * Vector
      *
      */
-    fn unit_complex_mul_vector(c: UnitComplex<f64>, v: Vector2<f64>, p: Point2<f64>) -> bool {
+    #[test]
+    fn unit_complex_mul_vector(c in unit_complex(), v in vector2(), p in point2()) {
         let r = c.to_rotation_matrix();
 
-        relative_eq!(c * v, r * v, epsilon = 1.0e-7) && relative_eq!(c * p, r * p, epsilon = 1.0e-7)
+        prop_assert!(relative_eq!(c * v, r * v, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(c * p, r * p, epsilon = 1.0e-7));
     }
 
     // Test that all operators (incl. all combinations of references) work.
     // See the top comment on `geometry/quaternion_ops.rs` for details on which operations are
     // supported.
+    #[test]
     #[cfg_attr(rustfmt, rustfmt_skip)]
     fn all_op_exist(
-        uc: UnitComplex<f64>,
-        v: Vector2<f64>,
-        p: Point2<f64>,
-        r: Rotation2<f64>
-    ) -> bool {
+        uc in unit_complex(),
+        v in vector2(),
+        p in point2(),
+        r in rotation2()
+    ) {
         let uv = Unit::new_normalize(v);
 
         let ucMuc = uc * uc;
@@ -112,7 +122,7 @@ quickcheck!(
         ucDr1 /= r;
         ucDr2 /= &r;
 
-        ucMuc1 == ucMuc
+        prop_assert!(ucMuc1 == ucMuc
             && ucMuc1 == ucMuc2
             && ucMr1 == ucMr
             && ucMr1 == ucMr2
@@ -146,6 +156,6 @@ quickcheck!(
             && ucMv == &uc * v
             && ucMuv == &uc * &uv
             && ucMuv == uc * &uv
-            && ucMuv == &uc * uv
+            && ucMuv == &uc * uv)
     }
 );

tests/lib.rs

@@ -12,9 +12,6 @@ extern crate approx;
 extern crate mint;
 extern crate nalgebra as na;
 extern crate num_traits as num;
-#[cfg(feature = "arbitrary")]
-#[macro_use]
-extern crate quickcheck;
 
 mod core;
 mod geometry;

tests/linalg/balancing.rs

@@ -1,26 +1,28 @@
-#![cfg(feature = "arbitrary")]
-
-use std::cmp;
+#![cfg(feature = "proptest-support")]
 
 use na::balancing;
-use na::{DMatrix, Matrix4};
+use na::DMatrix;
 
-quickcheck! {
-    fn balancing_parlett_reinsch(n: usize) -> bool {
-        let n = cmp::min(n, 10);
+use crate::proptest::*;
+use proptest::{prop_assert_eq, proptest};
+
+proptest! {
+    #[test]
+    fn balancing_parlett_reinsch(n in PROPTEST_MATRIX_DIM) {
         let m = DMatrix::<f64>::new_random(n, n);
         let mut balanced = m.clone();
         let d = balancing::balance_parlett_reinsch(&mut balanced);
         balancing::unbalance(&mut balanced, &d);
 
-        balanced == m
+        prop_assert_eq!(balanced, m);
     }
 
-    fn balancing_parlett_reinsch_static(m: Matrix4<f64>) -> bool {
+    #[test]
+    fn balancing_parlett_reinsch_static(m in matrix4()) {
         let mut balanced = m;
         let d = balancing::balance_parlett_reinsch(&mut balanced);
         balancing::unbalance(&mut balanced, &d);
 
-        balanced == m
+        prop_assert_eq!(balanced, m);
     }
 }

tests/linalg/bidiagonal.rs

@@ -1,63 +1,61 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 
 macro_rules! gen_tests(
-    ($module: ident, $scalar: ty) => {
+    ($module: ident, $scalar: expr) => {
         mod $module {
-            use na::{DMatrix, Matrix2, Matrix3x5, Matrix4, Matrix5x3};
             #[allow(unused_imports)]
             use crate::core::helper::{RandScalar, RandComplex};
 
-            quickcheck! {
-                fn bidiagonal(m: DMatrix<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
-                    if m.len() == 0  {
-                        return true;
-                    }
+            use crate::proptest::*;
+            use proptest::{prop_assert, proptest};
 
+            proptest! {
+                #[test]
+                fn bidiagonal(m in dmatrix_($scalar)) {
                     let bidiagonal = m.clone().bidiagonalize();
                     let (u, d, v_t) = bidiagonal.unpack();
 
-                    relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7))
                 }
 
-                fn bidiagonal_static_5_3(m: Matrix5x3<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn bidiagonal_static_5_3(m in matrix5x3_($scalar)) {
                     let bidiagonal = m.bidiagonalize();
                     let (u, d, v_t) = bidiagonal.unpack();
 
-                    relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7))
                 }
 
-                fn bidiagonal_static_3_5(m: Matrix3x5<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn bidiagonal_static_3_5(m in matrix3x5_($scalar)) {
                     let bidiagonal = m.bidiagonalize();
                     let (u, d, v_t) = bidiagonal.unpack();
 
-                    relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7))
                 }
 
-                fn bidiagonal_static_square(m: Matrix4<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn bidiagonal_static_square(m in matrix4_($scalar)) {
                     let bidiagonal = m.bidiagonalize();
                     let (u, d, v_t) = bidiagonal.unpack();
 
-                    relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7))
                 }
 
-                fn bidiagonal_static_square_2x2(m: Matrix2<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn bidiagonal_static_square_2x2(m in matrix2_($scalar)) {
                     let bidiagonal = m.bidiagonalize();
                     let (u, d, v_t) = bidiagonal.unpack();
 
-                    relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7))
                 }
             }
         }
     }
 );
 
-gen_tests!(complex, RandComplex<f64>);
-gen_tests!(f64, RandScalar<f64>);
+gen_tests!(complex, complex_f64());
+gen_tests!(f64, PROPTEST_F64);
 
 #[test]
 fn bidiagonal_identity() {

tests/linalg/cholesky.rs

@@ -1,4 +1,4 @@
-#![cfg(all(feature = "arbitrary", feature = "debug"))]
+#![cfg(all(feature = "proptest-support", feature = "debug"))]
 
 macro_rules! gen_tests(
     ($module: ident, $scalar: ty) => {
@@ -9,32 +9,30 @@ macro_rules! gen_tests(
             use rand::random;
             #[allow(unused_imports)]
             use crate::core::helper::{RandScalar, RandComplex};
-            use std::cmp;
 
-            quickcheck! {
-                fn cholesky(n: usize) -> bool {
-                    let m = RandomSDP::new(Dynamic::new(n.max(1).min(50)), || random::<$scalar>().0).unwrap();
+            use crate::proptest::*;
+            use proptest::{prop_assert, proptest};
+
+            proptest! {
+                #[test]
+                fn cholesky(n in PROPTEST_MATRIX_DIM) {
+                    let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
                     let l = m.clone().cholesky().unwrap().unpack();
-                    relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7));
                 }
 
-                fn cholesky_static(_m: RandomSDP<f64, U4>) -> bool {
+                #[test]
+                fn cholesky_static(_n in PROPTEST_MATRIX_DIM) {
                     let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
                     let chol = m.cholesky().unwrap();
                     let l    = chol.unpack();
 
-                    if !relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7) {
-                        false
-                    }
-                    else {
-                        true
-                    }
+                    prop_assert!(relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7));
                 }
 
-                fn cholesky_solve(n: usize, nb: usize) -> bool {
-                    let n = n.max(1).min(50);
+                #[test]
+                fn cholesky_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
                     let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
-                    let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
 
                     let chol = m.clone().cholesky().unwrap();
                     let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
@@ -43,11 +41,12 @@ macro_rules! gen_tests(
                     let sol1 = chol.solve(&b1);
                     let sol2 = chol.solve(&b2);
 
-                    relative_eq!(&m * &sol1, b1, epsilon = 1.0e-7) &&
-                    relative_eq!(&m * &sol2, b2, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(&m * &sol1, b1, epsilon = 1.0e-7));
+                    prop_assert!(relative_eq!(&m * &sol2, b2, epsilon = 1.0e-7));
                 }
 
-                fn cholesky_solve_static(_n: usize) -> bool {
+                #[test]
+                fn cholesky_solve_static(_n in PROPTEST_MATRIX_DIM) {
                     let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
                     let chol = m.clone().cholesky().unwrap();
                     let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
@@ -56,29 +55,32 @@ macro_rules! gen_tests(
                     let sol1 = chol.solve(&b1);
                     let sol2 = chol.solve(&b2);
 
-                    relative_eq!(m * sol1, b1, epsilon = 1.0e-7) &&
-                    relative_eq!(m * sol2, b2, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-7));
+                    prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-7));
                 }
 
-                fn cholesky_inverse(n: usize) -> bool {
-                    let m = RandomSDP::new(Dynamic::new(n.max(1).min(50)), || random::<$scalar>().0).unwrap();
+                #[test]
+                fn cholesky_inverse(n in PROPTEST_MATRIX_DIM) {
+                    let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
                     let m1 = m.clone().cholesky().unwrap().inverse();
                     let id1 = &m  * &m1;
                     let id2 = &m1 * &m;
 
-                    id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7)
+                    prop_assert!(id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7));
                 }
 
-                fn cholesky_inverse_static(_n: usize) -> bool {
+                #[test]
+                fn cholesky_inverse_static(_n in PROPTEST_MATRIX_DIM) {
                     let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
                     let m1 = m.clone().cholesky().unwrap().inverse();
                     let id1 = &m  * &m1;
                     let id2 = &m1 * &m;
 
-                    id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7)
+                    prop_assert!(id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7));
                 }
 
-                fn cholesky_rank_one_update(_n: usize) -> bool {
+                #[test]
+                fn cholesky_rank_one_update(_n in PROPTEST_MATRIX_DIM) {
                     let mut m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
                     let x = Vector4::<$scalar>::new_random().map(|e| e.0);
 
@@ -96,10 +98,11 @@ macro_rules! gen_tests(
                     // updates m manually
                     m.gerc(sigma_scalar, &x, &x, one); // m += sigma * x * x.adjoint()
 
-                    relative_eq!(m, m_chol_updated, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, m_chol_updated, epsilon = 1.0e-7));
                 }
 
-                fn cholesky_insert_column(n: usize) -> bool {
+                #[test]
+                fn cholesky_insert_column(n in PROPTEST_MATRIX_DIM) {
                     let n = n.max(1).min(10);
                     let j = random::<usize>() % n;
                     let m_updated = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
@@ -112,10 +115,11 @@ macro_rules! gen_tests(
                     let chol = m.clone().cholesky().unwrap().insert_column(j, col);
                     let m_chol_updated = chol.l() * chol.l().adjoint();
 
-                    relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7));
                 }
 
-                fn cholesky_remove_column(n: usize) -> bool {
+                #[test]
+                fn cholesky_remove_column(n in PROPTEST_MATRIX_DIM) {
                     let n = n.max(1).min(10);
                     let j = random::<usize>() % n;
                     let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
@@ -127,7 +131,7 @@ macro_rules! gen_tests(
                     // remove column from m
                     let m_updated = m.remove_column(j).remove_row(j);
 
-                    relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7));
                 }
             }
         }

tests/linalg/col_piv_qr.rs

@@ -22,133 +22,124 @@ fn col_piv_qr() {
     assert!(relative_eq!(m, qr, epsilon = 1.0e-7));
 }
 
-#[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
     macro_rules! gen_tests(
-        ($module: ident, $scalar: ty) => {
+        ($module: ident, $scalar: expr ,$scalar_type: ty) => {
             mod $module {
-                use na::{DMatrix, DVector, Matrix3x5, Matrix4, Matrix4x3, Matrix5x3, Vector4};
+                use na::{DMatrix, DVector, Matrix4x3, Vector4};
                 use std::cmp;
-                #[allow(unused_imports)]
-                use crate::core::helper::{RandScalar, RandComplex};
 
-                quickcheck! {
-                    fn col_piv_qr(m: DMatrix<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                #[allow(unused_imports)]
+                use crate::core::helper::{RandComplex, RandScalar};
+                use crate::proptest::*;
+                use proptest::{prop_assert, proptest};
+
+                proptest! {
+                    #[test]
+                    fn col_piv_qr(m in dmatrix_($scalar)) {
                         let col_piv_qr = m.clone().col_piv_qr();
                         let (q, r, p) = col_piv_qr.unpack();
                         let mut qr = &q * &r;
                         p.inv_permute_columns(&mut qr);
 
-                        println!("m: {}", m);
-                        println!("col_piv_qr: {}", &q * &r);
-
-                        relative_eq!(m, &qr, epsilon = 1.0e-7) &&
-                        q.is_orthogonal(1.0e-7)
+                        prop_assert!(relative_eq!(m, &qr, epsilon = 1.0e-7));
+                        prop_assert!(q.is_orthogonal(1.0e-7));
                     }
 
-                    fn col_piv_qr_static_5_3(m: Matrix5x3<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn col_piv_qr_static_5_3(m in matrix5x3_($scalar)) {
                         let col_piv_qr = m.col_piv_qr();
                         let (q, r, p) = col_piv_qr.unpack();
                         let mut qr = q * r;
                         p.inv_permute_columns(&mut qr);
 
-                        relative_eq!(m, qr, epsilon = 1.0e-7) &&
-                        q.is_orthogonal(1.0e-7)
+                        prop_assert!(relative_eq!(m, qr, epsilon = 1.0e-7));
+                        prop_assert!(q.is_orthogonal(1.0e-7));
                     }
 
-                    fn col_piv_qr_static_3_5(m: Matrix3x5<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn col_piv_qr_static_3_5(m in matrix3x5_($scalar)) {
                         let col_piv_qr = m.col_piv_qr();
                         let (q, r, p) = col_piv_qr.unpack();
                         let mut qr = q * r;
                         p.inv_permute_columns(&mut qr);
 
-                        relative_eq!(m, qr, epsilon = 1.0e-7) &&
-                        q.is_orthogonal(1.0e-7)
+                        prop_assert!(relative_eq!(m, qr, epsilon = 1.0e-7));
+                        prop_assert!(q.is_orthogonal(1.0e-7));
                     }
 
-                    fn col_piv_qr_static_square(m: Matrix4<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn col_piv_qr_static_square(m in matrix4_($scalar)) {
                         let col_piv_qr = m.col_piv_qr();
                         let (q, r, p) = col_piv_qr.unpack();
                         let mut qr = q * r;
                         p.inv_permute_columns(&mut qr);
 
-                        println!("{}{}{}{}", q, r, qr, m);
-
-                        relative_eq!(m, qr, epsilon = 1.0e-7) &&
-                        q.is_orthogonal(1.0e-7)
+                        prop_assert!(relative_eq!(m, qr, epsilon = 1.0e-7));
+                        prop_assert!(q.is_orthogonal(1.0e-7));
                     }
 
-                    fn col_piv_qr_solve(n: usize, nb: usize) -> bool {
+                    #[test]
+                    fn col_piv_qr_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
                         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::<$scalar>::new_random(n, n).map(|e| e.0);
+                            let m  = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
 
                             let col_piv_qr = m.clone().col_piv_qr();
-                            let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
-                            let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
+                            let b1 = DVector::<$scalar_type>::new_random(n).map(|e| e.0);
+                            let b2 = DMatrix::<$scalar_type>::new_random(n, nb).map(|e| e.0);
 
                             if col_piv_qr.is_invertible() {
                                 let sol1 = col_piv_qr.solve(&b1).unwrap();
                                 let sol2 = col_piv_qr.solve(&b2).unwrap();
 
-                                return relative_eq!(&m * sol1, b1, epsilon = 1.0e-6) &&
-                                    relative_eq!(&m * sol2, b2, epsilon = 1.0e-6)
+                                prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-6));
+                                prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-6));
+                            }
                         }
                     }
 
-                        return true;
-                    }
-
-                    fn col_piv_qr_solve_static(m: Matrix4<$scalar>) -> bool {
-                         let m = m.map(|e| e.0);
+                    #[test]
+                    fn col_piv_qr_solve_static(m in matrix4_($scalar)) {
                          let col_piv_qr = m.col_piv_qr();
-                         let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
-                         let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);
+                         let b1 = Vector4::<$scalar_type>::new_random().map(|e| e.0);
+                         let b2 = Matrix4x3::<$scalar_type>::new_random().map(|e| e.0);
 
                          if col_piv_qr.is_invertible() {
                              let sol1 = col_piv_qr.solve(&b1).unwrap();
                              let sol2 = col_piv_qr.solve(&b2).unwrap();
 
-                             relative_eq!(m * sol1, b1, epsilon = 1.0e-6) &&
-                             relative_eq!(m * sol2, b2, epsilon = 1.0e-6)
-                         }
-                         else {
-                             false
+                             prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-6));
+                             prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-6));
                          }
                     }
 
-                    fn col_piv_qr_inverse(n: usize) -> bool {
+                    #[test]
+                    fn col_piv_qr_inverse(n in PROPTEST_MATRIX_DIM) {
                         let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
-                        let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
+                        let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
 
                         if let Some(m1) = m.clone().col_piv_qr().try_inverse() {
                             let id1 = &m  * &m1;
                             let id2 = &m1 * &m;
 
-                            id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
-                        }
-                        else {
-                            true
+                            prop_assert!(id1.is_identity(1.0e-5));
+                            prop_assert!(id2.is_identity(1.0e-5));
                         }
                     }
 
-                    fn col_piv_qr_inverse_static(m: Matrix4<$scalar>) -> bool {
-                        let m  = m.map(|e| e.0);
+                    #[test]
+                    fn col_piv_qr_inverse_static(m in matrix4_($scalar)) {
                         let col_piv_qr = m.col_piv_qr();
 
                         if let Some(m1) = col_piv_qr.try_inverse() {
                             let id1 = &m  * &m1;
                             let id2 = &m1 * &m;
 
-                            id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
-                        }
-                        else {
-                            true
+                            prop_assert!(id1.is_identity(1.0e-5));
+                            prop_assert!(id2.is_identity(1.0e-5));
                         }
                     }
                 }
@@ -156,6 +147,6 @@ mod quickcheck_tests {
         }
     );
 
-    gen_tests!(complex, RandComplex<f64>);
-    gen_tests!(f64, RandScalar<f64>);
+    gen_tests!(complex, complex_f64(), RandComplex<f64>);
+    gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);
 }

tests/linalg/eigen.rs

@@ -1,66 +1,74 @@
 use na::DMatrix;
 
-#[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
     macro_rules! gen_tests(
-        ($module: ident, $scalar: ty) => {
+        ($module: ident, $scalar: expr, $scalar_type: ty) => {
             mod $module {
-                use na::{DMatrix, Matrix2, Matrix3, Matrix4};
+                use na::DMatrix;
                 #[allow(unused_imports)]
                 use crate::core::helper::{RandScalar, RandComplex};
                 use std::cmp;
 
-                quickcheck! {
-                    fn symmetric_eigen(n: usize) -> bool {
+                use crate::proptest::*;
+                use proptest::{prop_assert, proptest};
+
+                proptest! {
+                    #[test]
+                    fn symmetric_eigen(n in PROPTEST_MATRIX_DIM) {
                         let n      = cmp::max(1, cmp::min(n, 10));
-                        let m      = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0).hermitian_part();
+                        let m      = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0).hermitian_part();
                         let eig    = m.clone().symmetric_eigen();
                         let recomp = eig.recompose();
 
-                        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+                        prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5))
                     }
 
-                    fn symmetric_eigen_singular(n: usize) -> bool {
+                    #[test]
+                    fn symmetric_eigen_singular(n in PROPTEST_MATRIX_DIM) {
                         let n      = cmp::max(1, cmp::min(n, 10));
-                        let mut m  = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0).hermitian_part();
+                        let mut m  = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0).hermitian_part();
                         m.row_mut(n / 2).fill(na::zero());
                         m.column_mut(n / 2).fill(na::zero());
                         let eig    = m.clone().symmetric_eigen();
                         let recomp = eig.recompose();
 
-                        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+                        prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5))
                     }
 
-                    fn symmetric_eigen_static_square_4x4(m: Matrix4<$scalar>) -> bool {
-                        let m      = m.map(|e| e.0).hermitian_part();
+                    #[test]
+                    fn symmetric_eigen_static_square_4x4(m in matrix4_($scalar)) {
+                        let m      = m.hermitian_part();
                         let eig    = m.symmetric_eigen();
                         let recomp = eig.recompose();
 
-                        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+                        prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5))
                     }
 
-                    fn symmetric_eigen_static_square_3x3(m: Matrix3<$scalar>) -> bool {
-                        let m      = m.map(|e| e.0).hermitian_part();
+                    #[test]
+                    fn symmetric_eigen_static_square_3x3(m in matrix3_($scalar)) {
+                        let m      = m.hermitian_part();
                         let eig    = m.symmetric_eigen();
                         let recomp = eig.recompose();
 
-                        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+                        prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5))
                     }
 
-                    fn symmetric_eigen_static_square_2x2(m: Matrix2<$scalar>) -> bool {
-                        let m      = m.map(|e| e.0).hermitian_part();
+                    #[test]
+                    fn symmetric_eigen_static_square_2x2(m in matrix2_($scalar)) {
+                        let m      = m.hermitian_part();
                         let eig    = m.symmetric_eigen();
                         let recomp = eig.recompose();
 
-                        relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
+                        prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5))
                     }
                 }
             }
         }
     );
 
-    gen_tests!(complex, RandComplex<f64>);
-    gen_tests!(f64, RandScalar<f64>);
+    gen_tests!(complex, complex_f64(), RandComplex<f64>);
+    gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);
 }
 
 // Test proposed on the issue #176 of rulinalg.

tests/linalg/exp.rs

@@ -71,7 +71,6 @@ mod tests {
                 let m22 = ad_2.exp() * (delta * delta_2.cosh() + (d - a) * delta_2.sinh());
 
                 let f = Matrix2::new(m11, m12, m21, m22) / delta;
-                println!("a: {}", m);
                 assert!(relative_eq!(f, m.exp(), epsilon = 1.0e-7));
                 break;
             }

tests/linalg/full_piv_lu.rs

@@ -40,101 +40,96 @@ fn full_piv_lu_simple_with_pivot() {
 }
 
 #[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
+mod proptest_tests {
     macro_rules! gen_tests(
-    ($module: ident, $scalar: ty) => {
+    ($module: ident, $scalar: expr, $scalar_type: ty) => {
             mod $module {
                 use std::cmp;
                 use num::One;
-                use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
+                use na::{DMatrix, Matrix4x3, DVector, Vector4};
                 #[allow(unused_imports)]
                 use crate::core::helper::{RandScalar, RandComplex};
 
-                quickcheck! {
-                    fn full_piv_lu(m: DMatrix<$scalar>) -> bool {
-                        let mut m = m.map(|e| e.0);
-                        if m.len() == 0 {
-                            m = DMatrix::<$scalar>::new_random(1, 1).map(|e| e.0);
-                        }
+                use crate::proptest::*;
+                use proptest::{prop_assert, proptest};
 
+                proptest! {
+                    #[test]
+                    fn full_piv_lu(m in dmatrix_($scalar)) {
                         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)
+                        prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
                     }
 
-                    fn full_piv_lu_static_3_5(m: Matrix3x5<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn full_piv_lu_static_3_5(m in matrix3x5_($scalar)) {
                         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)
+                        prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
                     }
 
-                    fn full_piv_lu_static_5_3(m: Matrix5x3<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn full_piv_lu_static_5_3(m in matrix5x3_($scalar)) {
                         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)
+                        prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
                     }
 
-                    fn full_piv_lu_static_square(m: Matrix4<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn full_piv_lu_static_square(m in matrix4_($scalar)) {
                         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)
+                        prop_assert!(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::<$scalar>::new_random(n, n).map(|e| e.0);
+                    #[test]
+                    fn full_piv_lu_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
+                        let m  = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
 
                         let lu = m.clone().full_piv_lu();
-                            let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
-                            let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
+                        let b1 = DVector::<$scalar_type>::new_random(n).map(|e| e.0);
+                        let b2 = DMatrix::<$scalar_type>::new_random(n, nb).map(|e| e.0);
 
                         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))
+                        prop_assert!(sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6));
+                        prop_assert!(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6));
                     }
 
-                        return true;
-                    }
-
-                    fn full_piv_lu_solve_static(m: Matrix4<$scalar>) -> bool {
-                         let m = m.map(|e| e.0);
+                    #[test]
+                    fn full_piv_lu_solve_static(m in matrix4_($scalar)) {
                          let lu = m.full_piv_lu();
-                         let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
-                         let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);
+                         let b1 = Vector4::<$scalar_type>::new_random().map(|e| e.0);
+                         let b2 = Matrix4x3::<$scalar_type>::new_random().map(|e| e.0);
 
                          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))
+                         prop_assert!(sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6));
+                         prop_assert!(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6));
                     }
 
-                    fn full_piv_lu_inverse(n: usize) -> bool {
+                    #[test]
+                    fn full_piv_lu_inverse(n in PROPTEST_MATRIX_DIM) {
                         let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
-                        let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
+                        let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
 
                         let mut l = m.lower_triangle();
                         let mut u = m.upper_triangle();
@@ -148,21 +143,20 @@ mod quickcheck_tests {
                         let id1 = &m  * &m1;
                         let id2 = &m1 * &m;
 
-                        return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
+                        prop_assert!(id1.is_identity(1.0e-5));
+                        prop_assert!(id2.is_identity(1.0e-5));
                     }
 
-                    fn full_piv_lu_inverse_static(m: Matrix4<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn full_piv_lu_inverse_static(m in matrix4_($scalar)) {
                         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
+                            prop_assert!(id1.is_identity(1.0e-5));
+                            prop_assert!(id2.is_identity(1.0e-5));
                         }
                     }
                 }
@@ -170,8 +164,8 @@ mod quickcheck_tests {
         }
     );
 
-    gen_tests!(complex, RandComplex<f64>);
-    gen_tests!(f64, RandScalar<f64>);
+    gen_tests!(complex, complex_f64(), RandComplex<f64>);
+    gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);
 }
 
 /*

tests/linalg/hessenberg.rs

@@ -1,4 +1,4 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 
 use na::Matrix2;
 
@@ -11,40 +11,39 @@ fn hessenberg_simple() {
 }
 
 macro_rules! gen_tests(
-    ($module: ident, $scalar: ty) => {
+    ($module: ident, $scalar: expr) => {
          mod $module {
-            use na::{DMatrix, Matrix2, Matrix4};
-            use std::cmp;
             #[allow(unused_imports)]
             use crate::core::helper::{RandScalar, RandComplex};
 
-            quickcheck! {
-                fn hessenberg(n: usize) -> bool {
-                    let n = cmp::max(1, cmp::min(n, 50));
-                    let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
+            use crate::proptest::*;
+            use proptest::{prop_assert, proptest};
 
+            proptest! {
+                #[test]
+                fn hessenberg(m in dmatrix_($scalar)) {
                     let hess = m.clone().hessenberg();
                     let (p, h) = hess.unpack();
-                    relative_eq!(m, &p * h * p.adjoint(), epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, &p * h * p.adjoint(), epsilon = 1.0e-7))
                 }
 
-                fn hessenberg_static_mat2(m: Matrix2<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn hessenberg_static_mat2(m in matrix2_($scalar)) {
                     let hess = m.hessenberg();
                     let (p, h) = hess.unpack();
-                    relative_eq!(m, p * h * p.adjoint(), epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, p * h * p.adjoint(), epsilon = 1.0e-7))
                 }
 
-                fn hessenberg_static(m: Matrix4<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn hessenberg_static(m in matrix4_($scalar)) {
                     let hess = m.hessenberg();
                     let (p, h) = hess.unpack();
-                    relative_eq!(m, p * h * p.adjoint(), epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m, p * h * p.adjoint(), epsilon = 1.0e-7))
                 }
             }
          }
     }
 );
 
-gen_tests!(complex, RandComplex<f64>);
-gen_tests!(f64, RandScalar<f64>);
+gen_tests!(complex, complex_f64());
+gen_tests!(f64, PROPTEST_F64);

tests/linalg/lu.rs

@@ -38,103 +38,90 @@ fn lu_simple_with_pivot() {
     assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
 }
 
-#[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
-    #[allow(unused_imports)]
-    use crate::core::helper::{RandComplex, RandScalar};
-
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
     macro_rules! gen_tests(
-        ($module: ident, $scalar: ty) => {
+        ($module: ident, $scalar: expr, $scalar_type: ty) => {
             mod $module {
                 use std::cmp;
-                use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
+                use na::{DMatrix, Matrix4x3, DVector, Vector4};
                 #[allow(unused_imports)]
                 use crate::core::helper::{RandScalar, RandComplex};
+                use crate::proptest::*;
+                use proptest::{prop_assert, proptest};
 
-                quickcheck! {
-                    fn lu(m: DMatrix<$scalar>) -> bool {
-                        let mut m = m;
-                        if m.len() == 0 {
-                             m = DMatrix::<$scalar>::new_random(1, 1);
-                        }
-
-                        let m = m.map(|e| e.0);
-
+                proptest! {
+                    #[test]
+                    fn lu(m in dmatrix_($scalar)) {
                         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)
+                        prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
                     }
 
-                    fn lu_static_3_5(m: Matrix3x5<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn lu_static_3_5(m in matrix3x5_($scalar)) {
                         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)
+                        prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
                     }
 
-                    fn lu_static_5_3(m: Matrix5x3<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    fn lu_static_5_3(m in matrix5x3_($scalar)) {
                         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)
+                        prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
                     }
 
-                    fn lu_static_square(m: Matrix4<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn lu_static_square(m in matrix4_($scalar)) {
                         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)
+                        prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
                     }
 
-                    fn lu_solve(n: usize, nb: usize) -> bool {
-                        if n != 0 && nb != 0 {
+                    #[test]
+                    fn lu_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
                         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::<$scalar>::new_random(n, n).map(|e| e.0);
+                        let m  = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
 
                         let lu = m.clone().lu();
-                            let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
-                            let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
+                        let b1 = DVector::<$scalar_type>::new_random(n).map(|e| e.0);
+                        let b2 = DMatrix::<$scalar_type>::new_random(n, nb).map(|e| e.0);
 
                         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))
+                        prop_assert!(sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6));
+                        prop_assert!(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6));
                     }
 
-                        return true;
-                    }
-
-                    fn lu_solve_static(m: Matrix4<$scalar>) -> bool {
-                         let m = m.map(|e| e.0);
+                    #[test]
+                    fn lu_solve_static(m in matrix4_($scalar)) {
                          let lu = m.lu();
-                         let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
-                         let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);
+                         let b1 = Vector4::<$scalar_type>::new_random().map(|e| e.0);
+                         let b2 = Matrix4x3::<$scalar_type>::new_random().map(|e| e.0);
 
                          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))
+                         prop_assert!(sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6));
+                         prop_assert!(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::<$scalar>::new_random(n, n).map(|e| e.0);
-
+                    #[test]
+                    fn lu_inverse(m in dmatrix_($scalar)) {
                         let mut l = m.lower_triangle();
                         let mut u = m.upper_triangle();
 
@@ -147,21 +134,20 @@ mod quickcheck_tests {
                         let id1 = &m  * &m1;
                         let id2 = &m1 * &m;
 
-                        return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
+                        prop_assert!(id1.is_identity(1.0e-5));
+                        prop_assert!(id2.is_identity(1.0e-5));
                     }
 
-                    fn lu_inverse_static(m: Matrix4<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn lu_inverse_static(m in matrix4_($scalar)) {
                         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
+                            prop_assert!(id1.is_identity(1.0e-5));
+                            prop_assert!(id2.is_identity(1.0e-5));
                         }
                     }
                 }
@@ -169,6 +155,6 @@ mod quickcheck_tests {
         }
     );
 
-    gen_tests!(complex, RandComplex<f64>);
-    gen_tests!(f64, RandScalar<f64>);
+    gen_tests!(complex, complex_f64(), RandComplex<f64>);
+    gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);
 }

tests/linalg/qr.rs

@@ -1,126 +1,112 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 
 macro_rules! gen_tests(
-    ($module: ident, $scalar: ty) => {
+    ($module: ident, $scalar: expr, $scalar_type: ty) => {
         mod $module {
-            use na::{DMatrix, DVector, Matrix3x5, Matrix4, Matrix4x3, Matrix5x3, Vector4};
+            use na::{DMatrix, DVector, Matrix4x3, Vector4};
             use std::cmp;
             #[allow(unused_imports)]
             use crate::core::helper::{RandScalar, RandComplex};
+            use crate::proptest::*;
+            use proptest::{prop_assert, proptest};
 
-            quickcheck! {
-                fn qr(m: DMatrix<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+            proptest! {
+                #[test]
+                fn qr(m in dmatrix_($scalar)) {
                     let qr = m.clone().qr();
                     let q  = qr.q();
                     let r  = qr.r();
 
-                    println!("m: {}", m);
-                    println!("qr: {}", &q * &r);
-
-                    relative_eq!(m, &q * r, epsilon = 1.0e-7) &&
-                    q.is_orthogonal(1.0e-7)
+                    prop_assert!(relative_eq!(m, &q * r, epsilon = 1.0e-7));
+                    prop_assert!(q.is_orthogonal(1.0e-7));
                 }
 
-                fn qr_static_5_3(m: Matrix5x3<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn qr_static_5_3(m in matrix5x3_($scalar)) {
                     let qr = m.qr();
                     let q  = qr.q();
                     let r  = qr.r();
 
-                    relative_eq!(m, q * r, epsilon = 1.0e-7) &&
-                    q.is_orthogonal(1.0e-7)
+                    prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-7));
+                    prop_assert!(q.is_orthogonal(1.0e-7));
                 }
 
-                fn qr_static_3_5(m: Matrix3x5<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn qr_static_3_5(m in matrix3x5_($scalar)) {
                     let qr = m.qr();
                     let q  = qr.q();
                     let r  = qr.r();
 
-                    relative_eq!(m, q * r, epsilon = 1.0e-7) &&
-                    q.is_orthogonal(1.0e-7)
+                    prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-7));
+                    prop_assert!(q.is_orthogonal(1.0e-7));
                 }
 
-                fn qr_static_square(m: Matrix4<$scalar>) -> bool {
-                    let m = m.map(|e| e.0);
+                #[test]
+                fn qr_static_square(m in matrix4_($scalar)) {
                     let qr = m.qr();
                     let q  = qr.q();
                     let r  = qr.r();
 
-                    println!("{}{}{}{}", q, r, q * r, m);
-
-                    relative_eq!(m, q * r, epsilon = 1.0e-7) &&
-                    q.is_orthogonal(1.0e-7)
+                    prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-7));
+                    prop_assert!(q.is_orthogonal(1.0e-7));
                 }
 
-                fn qr_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::<$scalar>::new_random(n, n).map(|e| e.0);
+                #[test]
+                fn qr_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
+                    let m  = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
 
                     let qr = m.clone().qr();
-                        let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
-                        let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
+                    let b1 = DVector::<$scalar_type>::new_random(n).map(|e| e.0);
+                    let b2 = DMatrix::<$scalar_type>::new_random(n, nb).map(|e| e.0);
 
                     if qr.is_invertible() {
                         let sol1 = qr.solve(&b1).unwrap();
                         let sol2 = qr.solve(&b2).unwrap();
 
-                            return relative_eq!(&m * sol1, b1, epsilon = 1.0e-6) &&
-                                relative_eq!(&m * sol2, b2, epsilon = 1.0e-6)
+                        prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-6));
+                        prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-6));
                     }
                 }
 
-                    return true;
-                }
-
-                fn qr_solve_static(m: Matrix4<$scalar>) -> bool {
-                     let m = m.map(|e| e.0);
+                #[test]
+                fn qr_solve_static(m in matrix4_($scalar)) {
                      let qr = m.qr();
-                     let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
-                     let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);
+                     let b1 = Vector4::<$scalar_type>::new_random().map(|e| e.0);
+                     let b2 = Matrix4x3::<$scalar_type>::new_random().map(|e| e.0);
 
                      if qr.is_invertible() {
                          let sol1 = qr.solve(&b1).unwrap();
                          let sol2 = qr.solve(&b2).unwrap();
 
-                         relative_eq!(m * sol1, b1, epsilon = 1.0e-6) &&
-                         relative_eq!(m * sol2, b2, epsilon = 1.0e-6)
-                     }
-                     else {
-                         false
+                         prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-6));
+                         prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-6));
                      }
                 }
 
-                fn qr_inverse(n: usize) -> bool {
+                #[test]
+                fn qr_inverse(n in PROPTEST_MATRIX_DIM) {
                     let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
-                    let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
+                    let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
 
                     if let Some(m1) = m.clone().qr().try_inverse() {
                         let id1 = &m  * &m1;
                         let id2 = &m1 * &m;
 
-                        id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
-                    }
-                    else {
-                        true
+                        prop_assert!(id1.is_identity(1.0e-5));
+                        prop_assert!(id2.is_identity(1.0e-5));
                     }
                 }
 
-                fn qr_inverse_static(m: Matrix4<$scalar>) -> bool {
-                    let m  = m.map(|e| e.0);
+                #[test]
+                fn qr_inverse_static(m in matrix4_($scalar)) {
                     let qr = m.qr();
 
                     if let Some(m1) = qr.try_inverse() {
                         let id1 = &m  * &m1;
                         let id2 = &m1 * &m;
 
-                        id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
-                    }
-                    else {
-                        true
+                        prop_assert!(id1.is_identity(1.0e-5));
+                        prop_assert!(id2.is_identity(1.0e-5));
                     }
                 }
             }
@@ -128,5 +114,5 @@ macro_rules! gen_tests(
     }
 );
 
-gen_tests!(complex, RandComplex<f64>);
-gen_tests!(f64, RandScalar<f64>);
+gen_tests!(complex, complex_f64(), RandComplex<f64>);
+gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);

tests/linalg/schur.rs

@@ -13,72 +13,47 @@ fn schur_simpl_mat3() {
     assert!(relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7));
 }
 
-#[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
     macro_rules! gen_tests(
-        ($module: ident, $scalar: ty) => {
+        ($module: ident, $scalar: expr) => {
             mod $module {
-                use std::cmp;
-                use na::{DMatrix, Matrix2, Matrix3, Matrix4};
                 #[allow(unused_imports)]
                 use crate::core::helper::{RandScalar, RandComplex};
+                use crate::proptest::*;
+                use proptest::{prop_assert, proptest};
 
-                quickcheck! {
-                    fn schur(n: usize) -> bool {
-                        let n = cmp::max(1, cmp::min(n, 10));
-                        let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
-
+                proptest! {
+                    #[test]
+                    fn schur(m in dmatrix_($scalar)) {
                         let (vecs, vals) = m.clone().schur().unpack();
-
-                        if !relative_eq!(&vecs * &vals * vecs.adjoint(), m, epsilon = 1.0e-7) {
-                            println!("{:.5}{:.5}", m, &vecs * &vals * vecs.adjoint());
+                        prop_assert!(relative_eq!(&vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7));
                     }
 
-                        relative_eq!(&vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7)
-                    }
-
-                    fn schur_static_mat2(m: Matrix2<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn schur_static_mat2(m in matrix2_($scalar)) {
                         let (vecs, vals) = m.clone().schur().unpack();
-
-                        let ok = relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7);
-                        if !ok {
-                            println!("Vecs: {:.5} Vals: {:.5}", vecs, vals);
-                            println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.adjoint());
-                        }
-                        ok
+                        prop_assert!(relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7));
                     }
 
-                    fn schur_static_mat3(m: Matrix3<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn schur_static_mat3(m in matrix3_($scalar)) {
                         let (vecs, vals) = m.clone().schur().unpack();
-
-                        let ok = relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7);
-                        if !ok {
-                            println!("Vecs: {:.5} Vals: {:.5}", vecs, vals);
-                            println!("{:.5}{:.5}", m, &vecs * &vals * vecs.adjoint());
-                        }
-                        ok
+                        prop_assert!(relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7));
                     }
 
-                    fn schur_static_mat4(m: Matrix4<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn schur_static_mat4(m in matrix4_($scalar)) {
                         let (vecs, vals) = m.clone().schur().unpack();
-
-                        let ok = relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7);
-                        if !ok {
-                            println!("{:.5}{:.5}", m, &vecs * &vals * vecs.adjoint());
-                        }
-
-                        ok
+                        prop_assert!(relative_eq!(vecs * vals * vecs.adjoint(), m, epsilon = 1.0e-7));
                     }
                 }
             }
         }
     );
 
-    gen_tests!(complex, RandComplex<f64>);
-    gen_tests!(f64, RandScalar<f64>);
+    gen_tests!(complex, complex_f64());
+    gen_tests!(f64, PROPTEST_F64);
 }
 
 #[test]

tests/linalg/solve.rs

@@ -1,11 +1,13 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 
 macro_rules! gen_tests(
-    ($module: ident, $scalar: ty) => {
+    ($module: ident, $scalar: expr) => {
         mod $module {
-            use na::{Matrix4, Matrix4x5, ComplexField};
+            use na::{Matrix4, ComplexField};
             #[allow(unused_imports)]
             use crate::core::helper::{RandScalar, RandComplex};
+            use crate::proptest::*;
+            use proptest::{prop_assert, proptest};
 
             fn unzero_diagonal<N: ComplexField>(a: &mut Matrix4<N>) {
                 for i in 0..4 {
@@ -15,50 +17,50 @@ macro_rules! gen_tests(
                 }
             }
 
-            quickcheck! {
-                fn solve_lower_triangular(a: Matrix4<$scalar>, b: Matrix4x5<$scalar>) -> bool {
-                    let b = b.map(|e| e.0);
-                    let mut a = a.map(|e| e.0);
+            proptest! {
+                #[test]
+                fn solve_lower_triangular(a in matrix4_($scalar), b in matrix4x5_($scalar)) {
+                    let mut a = a;
                     unzero_diagonal(&mut a);
                     let tri = a.lower_triangle();
                     let x   = a.solve_lower_triangular(&b).unwrap();
 
-                    relative_eq!(tri * x, b, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(tri * x, b, epsilon = 1.0e-7))
                 }
 
-                fn solve_upper_triangular(a: Matrix4<$scalar>, b: Matrix4x5<$scalar>) -> bool {
-                    let b = b.map(|e| e.0);
-                    let mut a = a.map(|e| e.0);
+                #[test]
+                fn solve_upper_triangular(a in matrix4_($scalar), b in matrix4x5_($scalar)) {
+                    let mut a = a;
                     unzero_diagonal(&mut a);
                     let tri = a.upper_triangle();
                     let x   = a.solve_upper_triangular(&b).unwrap();
 
-                    relative_eq!(tri * x, b, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(tri * x, b, epsilon = 1.0e-7))
                 }
 
-                fn tr_solve_lower_triangular(a: Matrix4<$scalar>, b: Matrix4x5<$scalar>) -> bool {
-                    let b = b.map(|e| e.0);
-                    let mut a = a.map(|e| e.0);
+                #[test]
+                fn tr_solve_lower_triangular(a in matrix4_($scalar), b in matrix4x5_($scalar)) {
+                    let mut a = a;
                     unzero_diagonal(&mut a);
                     let tri = a.lower_triangle();
                     let x   = a.tr_solve_lower_triangular(&b).unwrap();
 
-                    relative_eq!(tri.transpose() * x, b, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(tri.transpose() * x, b, epsilon = 1.0e-7))
                 }
 
-                fn tr_solve_upper_triangular(a: Matrix4<$scalar>, b: Matrix4x5<$scalar>) -> bool {
-                    let b = b.map(|e| e.0);
-                    let mut a = a.map(|e| e.0);
+                #[test]
+                fn tr_solve_upper_triangular(a in matrix4_($scalar), b in matrix4x5_($scalar)) {
+                    let mut a = a;
                     unzero_diagonal(&mut a);
                     let tri = a.upper_triangle();
                     let x   = a.tr_solve_upper_triangular(&b).unwrap();
 
-                    relative_eq!(tri.transpose() * x, b, epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(tri.transpose() * x, b, epsilon = 1.0e-7))
                 }
             }
         }
     }
 );
 
-gen_tests!(complex, RandComplex<f64>);
-gen_tests!(f64, RandScalar<f64>);
+gen_tests!(complex, complex_f64());
+gen_tests!(f64, PROPTEST_F64);

tests/linalg/svd.rs

@@ -1,162 +1,143 @@
 use na::{DMatrix, Matrix6};
 
-#[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
     macro_rules! gen_tests(
-        ($module: ident, $scalar: ty) => {
+        ($module: ident, $scalar: expr, $scalar_type: ty) => {
             mod $module {
                 use na::{
-                    DMatrix, DVector, Matrix2, Matrix2x5, Matrix3, Matrix3x5, Matrix4, Matrix5x2, Matrix5x3,
+                    DMatrix, DVector, Matrix2, Matrix3, Matrix4,
                     ComplexField
                 };
                 use std::cmp;
                 #[allow(unused_imports)]
                 use crate::core::helper::{RandScalar, RandComplex};
+                use crate::proptest::*;
+                use proptest::{prop_assert, proptest};
 
-                quickcheck! {
-                    fn svd(m: DMatrix<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
-                        if m.len() > 0 {
+                proptest! {
+                    #[test]
+                    fn svd(m in dmatrix_($scalar)) {
                         let svd = m.clone().svd(true, true);
                         let recomp_m = svd.clone().recompose().unwrap();
                         let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
                         let ds = DMatrix::from_diagonal(&s.map(|e| ComplexField::from_real(e)));
 
-                            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
-                        }
+                        prop_assert!(s.iter().all(|e| *e >= 0.0));
+                        prop_assert!(relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5));
+                        prop_assert!(relative_eq!(m, recomp_m, epsilon = 1.0e-5));
                     }
 
-                    fn svd_static_5_3(m: Matrix5x3<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn svd_static_5_3(m in matrix5x3_($scalar)) {
                         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.map(|e| ComplexField::from_real(e)));
 
-                        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)
+                        prop_assert!(s.iter().all(|e| *e >= 0.0));
+                        prop_assert!(relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5));
+                        prop_assert!(u.is_orthogonal(1.0e-5));
+                        prop_assert!(v_t.is_orthogonal(1.0e-5));
                     }
 
-                    fn svd_static_5_2(m: Matrix5x2<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn svd_static_5_2(m in matrix5x2_($scalar)) {
                         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.map(|e| ComplexField::from_real(e)));
 
-                        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)
+                        prop_assert!(s.iter().all(|e| *e >= 0.0));
+                        prop_assert!(relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5));
+                        prop_assert!(u.is_orthogonal(1.0e-5));
+                        prop_assert!(v_t.is_orthogonal(1.0e-5));
                     }
 
-                    fn svd_static_3_5(m: Matrix3x5<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn svd_static_3_5(m in matrix3x5_($scalar)) {
                         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.map(|e| ComplexField::from_real(e)));
 
-                        s.iter().all(|e| *e >= 0.0) &&
-                        relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
+                        prop_assert!(s.iter().all(|e| *e >= 0.0));
+                        prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5));
                     }
 
-                    fn svd_static_2_5(m: Matrix2x5<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn svd_static_2_5(m in matrix2x5_($scalar)) {
                         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.map(|e| ComplexField::from_real(e)));
 
-                        s.iter().all(|e| *e >= 0.0) &&
-                        relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
+                        prop_assert!(s.iter().all(|e| *e >= 0.0));
+                        prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5));
                     }
 
-                    fn svd_static_square(m: Matrix4<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn svd_static_square(m in matrix4_($scalar)) {
                         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.map(|e| ComplexField::from_real(e)));
 
-                        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)
+                        prop_assert!(s.iter().all(|e| *e >= 0.0));
+                        prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5));
+                        prop_assert!(u.is_orthogonal(1.0e-5));
+                        prop_assert!(v_t.is_orthogonal(1.0e-5));
                     }
 
-                    fn svd_static_square_2x2(m: Matrix2<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
+                    #[test]
+                    fn svd_static_square_2x2(m in matrix2_($scalar)) {
                         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.map(|e| ComplexField::from_real(e)));
 
-                        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)
+                        prop_assert!(s.iter().all(|e| *e >= 0.0));
+                        prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5));
+                        prop_assert!(u.is_orthogonal(1.0e-5));
+                        prop_assert!(v_t.is_orthogonal(1.0e-5));
                     }
 
-                    fn svd_pseudo_inverse(m: DMatrix<$scalar>) -> bool {
-                        let m = m.map(|e| e.0);
-
-                        if m.len() > 0 {
+                    #[test]
+                    fn svd_pseudo_inverse(m in dmatrix_($scalar)) {
                         let svd = m.clone().svd(true, true);
                         let pinv = svd.pseudo_inverse(1.0e-10).unwrap();
 
                         if m.nrows() > m.ncols() {
-                                (pinv * m).is_identity(1.0e-5)
-                            }
-                            else {
-                                (m * pinv).is_identity(1.0e-5)
-                            }
-                        }
-                        else {
-                            true
+                            prop_assert!((pinv * m).is_identity(1.0e-5))
+                        } else {
+                            prop_assert!((m * pinv).is_identity(1.0e-5))
                         }
                     }
 
-                    fn svd_solve(n: usize, nb: usize) -> bool {
+                    #[test]
+                    fn svd_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
                         let n = cmp::max(1, cmp::min(n, 10));
                         let nb = cmp::min(nb, 10);
-                        let m  = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
+                        let m  = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
 
                         let svd = m.clone().svd(true, true);
 
                         if svd.rank(1.0e-7) == n {
-                            let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
-                            let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
+                            let b1 = DVector::<$scalar_type>::new_random(n).map(|e| e.0);
+                            let b2 = DMatrix::<$scalar_type>::new_random(n, nb).map(|e| e.0);
 
                             let sol1 = svd.solve(&b1, 1.0e-7).unwrap();
                             let sol2 = svd.solve(&b2, 1.0e-7).unwrap();
 
                             let recomp = svd.recompose().unwrap();
-                            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;
+                            prop_assert!(relative_eq!(m, recomp, epsilon = 1.0e-6));
+                            prop_assert!(relative_eq!(&m * &sol1, b1, epsilon = 1.0e-6));
+                            prop_assert!(relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6));
                         }
-                            if !relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6) {
-                                println!("Problem 2: {:.6}{:.6}", b2, &m * sol2);
-                                return false;
-                            }
-                        }
-
-                        true
                     }
                 }
             }
         }
     );
 
-    gen_tests!(complex, RandComplex<f64>);
-    gen_tests!(f64, RandScalar<f64>);
+    gen_tests!(complex, complex_f64(), RandComplex<f64>);
+    gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);
 }
 
 // Test proposed on the issue #176 of rulinalg.
@@ -303,31 +284,31 @@ fn svd_identity() {
 #[rustfmt::skip]
 fn svd_with_delimited_subproblem() {
     let mut m = DMatrix::<f64>::from_element(10, 10, 0.0);
-    m[(0,0)] = 1.0;  m[(0,1)] = 2.0;
-    m[(1,1)] = 0.0;  m[(1,2)] = 3.0;
-    m[(2,2)] = 4.0;  m[(2,3)] = 5.0;
-    m[(3,3)] = 6.0;  m[(3,4)] = 0.0;
-    m[(4,4)] = 8.0;  m[(3,5)] = 9.0;
-    m[(5,5)] = 10.0; m[(3,6)] = 11.0;
-    m[(6,6)] = 12.0; m[(3,7)] = 12.0;
-    m[(7,7)] = 14.0; m[(3,8)] = 13.0;
-    m[(8,8)] = 16.0; m[(3,9)] = 17.0;
-    m[(9,9)] = 18.0;
+    m[(0, 0)] = 1.0;  m[(0, 1)] = 2.0;
+    m[(1, 1)] = 0.0;  m[(1, 2)] = 3.0;
+    m[(2, 2)] = 4.0;  m[(2, 3)] = 5.0;
+    m[(3, 3)] = 6.0;  m[(3, 4)] = 0.0;
+    m[(4, 4)] = 8.0;  m[(3, 5)] = 9.0;
+    m[(5, 5)] = 10.0; m[(3, 6)] = 11.0;
+    m[(6, 6)] = 12.0; m[(3, 7)] = 12.0;
+    m[(7, 7)] = 14.0; m[(3, 8)] = 13.0;
+    m[(8, 8)] = 16.0; m[(3, 9)] = 17.0;
+    m[(9, 9)] = 18.0;
     let svd = m.clone().svd(true, true);
     assert_relative_eq!(m, svd.recompose().unwrap(), epsilon = 1.0e-7);
 
     // Rectangular versions.
     let mut m = DMatrix::<f64>::from_element(15, 10, 0.0);
-    m[(0,0)] = 1.0;  m[(0,1)] = 2.0;
-    m[(1,1)] = 0.0;  m[(1,2)] = 3.0;
-    m[(2,2)] = 4.0;  m[(2,3)] = 5.0;
-    m[(3,3)] = 6.0;  m[(3,4)] = 0.0;
-    m[(4,4)] = 8.0;  m[(3,5)] = 9.0;
-    m[(5,5)] = 10.0; m[(3,6)] = 11.0;
-    m[(6,6)] = 12.0; m[(3,7)] = 12.0;
-    m[(7,7)] = 14.0; m[(3,8)] = 13.0;
-    m[(8,8)] = 16.0; m[(3,9)] = 17.0;
-    m[(9,9)] = 18.0;
+    m[(0, 0)] = 1.0;  m[(0, 1)] = 2.0;
+    m[(1, 1)] = 0.0;  m[(1, 2)] = 3.0;
+    m[(2, 2)] = 4.0;  m[(2, 3)] = 5.0;
+    m[(3, 3)] = 6.0;  m[(3, 4)] = 0.0;
+    m[(4, 4)] = 8.0;  m[(3, 5)] = 9.0;
+    m[(5, 5)] = 10.0; m[(3, 6)] = 11.0;
+    m[(6, 6)] = 12.0; m[(3, 7)] = 12.0;
+    m[(7, 7)] = 14.0; m[(3, 8)] = 13.0;
+    m[(8, 8)] = 16.0; m[(3, 9)] = 17.0;
+    m[(9, 9)] = 18.0;
     let svd = m.clone().svd(true, true);
     assert_relative_eq!(m, svd.recompose().unwrap(), epsilon = 1.0e-7);
 

tests/linalg/tridiagonal.rs

@@ -1,54 +1,56 @@
-#![cfg(feature = "arbitrary")]
+#![cfg(feature = "proptest-support")]
 
 macro_rules! gen_tests(
-    ($module: ident, $scalar: ty) => {
+    ($module: ident, $scalar: expr) => {
             mod $module {
-            use std::cmp;
-
-            use na::{DMatrix, Matrix2, Matrix4};
             #[allow(unused_imports)]
             use crate::core::helper::{RandScalar, RandComplex};
+            use crate::proptest::*;
+            use proptest::{prop_assert, proptest};
 
-            quickcheck! {
-                fn symm_tridiagonal(n: usize) -> bool {
-                    let n = cmp::max(1, cmp::min(n, 50));
-                    let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0).hermitian_part();
+            proptest! {
+                #[test]
+                fn symm_tridiagonal(m in dmatrix_($scalar)) {
+                    let m = &m * m.adjoint();
                     let tri = m.clone().symmetric_tridiagonalize();
                     let recomp = tri.recompose();
 
-                    relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7));
                 }
 
-                fn symm_tridiagonal_singular(n: usize) -> bool {
-                    let n = cmp::max(1, cmp::min(n, 4));
-                    let mut m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0).hermitian_part();
+                #[test]
+                fn symm_tridiagonal_singular(m in dmatrix_($scalar)) {
+                    let mut m = &m * m.adjoint();
+                    let n = m.nrows();
                     m.row_mut(n / 2).fill(na::zero());
                     m.column_mut(n / 2).fill(na::zero());
                     let tri = m.clone().symmetric_tridiagonalize();
                     let recomp = tri.recompose();
 
-                    relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7));
                 }
 
-                fn symm_tridiagonal_static_square(m: Matrix4<$scalar>) -> bool {
-                    let m = m.map(|e| e.0).hermitian_part();
+                #[test]
+                fn symm_tridiagonal_static_square(m in matrix4_($scalar)) {
+                    let m = m.hermitian_part();
                     let tri = m.symmetric_tridiagonalize();
                     let recomp = tri.recompose();
 
-                    relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7));
                 }
 
-                fn symm_tridiagonal_static_square_2x2(m: Matrix2<$scalar>) -> bool {
-                    let m = m.map(|e| e.0).hermitian_part();
+                #[test]
+                fn symm_tridiagonal_static_square_2x2(m in matrix2_($scalar)) {
+                    let m = m.hermitian_part();
                     let tri = m.symmetric_tridiagonalize();
                     let recomp = tri.recompose();
 
-                    relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
+                    prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7));
                 }
             }
         }
     }
 );
 
-gen_tests!(complex, RandComplex<f64>);
-gen_tests!(f64, RandScalar<f64>);
+gen_tests!(complex, complex_f64());
+gen_tests!(f64, PROPTEST_F64);

tests/linalg/udu.rs

@@ -8,7 +8,7 @@ fn udu_simple() {
        -1.0,  2.0, -1.0,
         0.0, -1.0,  2.0);
 
-    let udu = m.udu();
+    let udu = m.udu().unwrap();
     
     // Rebuild
     let p = udu.u * udu.d_matrix() * udu.u.transpose();
@@ -23,50 +23,54 @@ fn udu_non_sym_panic() {
     let m = Matrix3::new(
         2.0, -1.0,  0.0,
         1.0, -2.0,  3.0,
-       -2.0,  1.0,  0.0);
+       -2.0,  1.0,  0.3);
 
-    let udu = m.udu();
+    let udu = m.udu().unwrap();
     // Rebuild
     let p = udu.u * udu.d_matrix() * udu.u.transpose();
 
     assert!(relative_eq!(m, p, epsilon = 3.0e-16));
 }
 
-#[cfg(feature = "arbitrary")]
-mod quickcheck_tests {
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
     #[allow(unused_imports)]
     use crate::core::helper::{RandComplex, RandScalar};
 
     macro_rules! gen_tests(
-        ($module: ident, $scalar: ty) => {
+        ($module: ident, $scalar: expr) => {
             mod $module {
-                use na::{DMatrix, Matrix4};
                 #[allow(unused_imports)]
                 use crate::core::helper::{RandScalar, RandComplex};
+                use crate::proptest::*;
+                use proptest::{prop_assert, proptest};
 
-                quickcheck! {
-                    fn udu(n: usize) -> bool {
-                        let n = std::cmp::max(1, std::cmp::min(n, 10));
-                        let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0).hermitian_part();
+                proptest! {
+                    #[test]
+                    fn udu(m in dmatrix_($scalar)) {
+                        let m = &m * m.adjoint();
 
-                        let udu = m.clone().udu();
+                        if let Some(udu) = m.clone().udu() {
                             let p = &udu.u * &udu.d_matrix() * &udu.u.transpose();
+                            println!("m: {}, p: {}", m, p);
 
-                        relative_eq!(m, p, epsilon = 1.0e-7)
+                            prop_assert!(relative_eq!(m, p, epsilon = 1.0e-7));
+                        }
                     }
 
-                    fn udu_static(m: Matrix4<$scalar>) -> bool {
-                        let m = m.map(|e| e.0).hermitian_part();
+                    #[test]
+                    fn udu_static(m in matrix4_($scalar)) {
+                        let m = m.hermitian_part();
 
-                        let udu = m.udu();
+                        if let Some(udu) = m.udu() {
                             let p = udu.u * udu.d_matrix() * udu.u.transpose();
-
-                        relative_eq!(m, p, epsilon = 1.0e-7)
+                            prop_assert!(relative_eq!(m, p, epsilon = 1.0e-7));
+                        }
                     }
                 }
             }
         }
     );
 
-    gen_tests!(f64, RandScalar<f64>);
+    gen_tests!(f64, PROPTEST_F64);
 }

tests/proptest/mod.rs

@@ -1,10 +1,161 @@
 //! Tests for proptest-related functionality.
+use nalgebra::allocator::Allocator;
 use nalgebra::base::dimension::*;
-use nalgebra::proptest::{matrix, DimRange, MatrixStrategy};
-use nalgebra::{DMatrix, DVector, Dim, Matrix3, MatrixMN, Vector3};
+use nalgebra::proptest::{DimRange, MatrixStrategy};
+use nalgebra::{
+    DMatrix, DVector, DefaultAllocator, Dim, DualQuaternion, Isometry2, Isometry3, Matrix3,
+    MatrixMN, Point2, Point3, Quaternion, Rotation2, Rotation3, Scalar, Similarity3, Translation2,
+    Translation3, UnitComplex, UnitDualQuaternion, UnitQuaternion, Vector3, U2, U3, U4, U7, U8,
+};
+use num_complex::Complex;
 use proptest::prelude::*;
-use proptest::strategy::ValueTree;
+use proptest::strategy::{Strategy, ValueTree};
 use proptest::test_runner::TestRunner;
+use std::ops::RangeInclusive;
+
+pub const PROPTEST_MATRIX_DIM: RangeInclusive<usize> = 1..=20;
+pub const PROPTEST_F64: RangeInclusive<f64> = -100.0..=100.0;
+
+pub use nalgebra::proptest::{matrix, vector};
+
+pub fn point2() -> impl Strategy<Value = Point2<f64>> {
+    vector2().prop_map(|v| Point2::from(v))
+}
+
+pub fn point3() -> impl Strategy<Value = Point3<f64>> {
+    vector3().prop_map(|v| Point3::from(v))
+}
+
+pub fn translation2() -> impl Strategy<Value = Translation2<f64>> {
+    vector2().prop_map(|v| Translation2::from(v))
+}
+
+pub fn translation3() -> impl Strategy<Value = Translation3<f64>> {
+    vector3().prop_map(|v| Translation3::from(v))
+}
+
+pub fn rotation2() -> impl Strategy<Value = Rotation2<f64>> {
+    PROPTEST_F64.prop_map(|v| Rotation2::new(v))
+}
+
+pub fn rotation3() -> impl Strategy<Value = Rotation3<f64>> {
+    vector3().prop_map(|v| Rotation3::new(v))
+}
+
+pub fn unit_complex() -> impl Strategy<Value = UnitComplex<f64>> {
+    PROPTEST_F64.prop_map(|v| UnitComplex::new(v))
+}
+
+pub fn isometry2() -> impl Strategy<Value = Isometry2<f64>> {
+    vector3().prop_map(|v| Isometry2::new(v.xy(), v.z))
+}
+
+pub fn isometry3() -> impl Strategy<Value = Isometry3<f64>> {
+    vector6().prop_map(|v| Isometry3::new(v.xyz(), Vector3::new(v.w, v.a, v.b)))
+}
+
+// pub fn similarity2() -> impl Strategy<Value = Similarity2<f64>> {
+//     vector4().prop_map(|v| Similarity2::new(v.xy(), v.z, v.w))
+// }
+
+pub fn similarity3() -> impl Strategy<Value = Similarity3<f64>> {
+    vector(PROPTEST_F64, U7)
+        .prop_map(|v| Similarity3::new(v.xyz(), Vector3::new(v[3], v[4], v[5]), v[6]))
+}
+
+pub fn unit_dual_quaternion() -> impl Strategy<Value = UnitDualQuaternion<f64>> {
+    isometry3().prop_map(|iso| UnitDualQuaternion::from_isometry(&iso))
+}
+
+pub fn dual_quaternion() -> impl Strategy<Value = DualQuaternion<f64>> {
+    vector(PROPTEST_F64, U8).prop_map(|v| {
+        DualQuaternion::from_real_and_dual(
+            Quaternion::new(v[0], v[1], v[2], v[3]),
+            Quaternion::new(v[4], v[5], v[6], v[7]),
+        )
+    })
+}
+
+pub fn quaternion() -> impl Strategy<Value = Quaternion<f64>> {
+    vector4().prop_map(|v| Quaternion::from(v))
+}
+
+pub fn unit_quaternion() -> impl Strategy<Value = UnitQuaternion<f64>> {
+    vector3().prop_map(|v| UnitQuaternion::new(v))
+}
+
+pub fn complex_f64() -> impl Strategy<Value = Complex<f64>> + Clone {
+    vector(PROPTEST_F64, U2).prop_map(|v| Complex::new(v.x, v.y))
+}
+
+pub fn dmatrix() -> impl Strategy<Value = DMatrix<f64>> {
+    matrix(PROPTEST_F64, PROPTEST_MATRIX_DIM, PROPTEST_MATRIX_DIM)
+}
+
+pub fn dvector() -> impl Strategy<Value = DVector<f64>> {
+    vector(PROPTEST_F64, PROPTEST_MATRIX_DIM)
+}
+
+pub fn dmatrix_<ScalarStrategy>(
+    scalar_strategy: ScalarStrategy,
+) -> impl Strategy<Value = DMatrix<ScalarStrategy::Value>>
+where
+    ScalarStrategy: Strategy + Clone + 'static,
+    ScalarStrategy::Value: Scalar,
+    DefaultAllocator: Allocator<ScalarStrategy::Value, Dynamic, Dynamic>,
+{
+    matrix(scalar_strategy, PROPTEST_MATRIX_DIM, PROPTEST_MATRIX_DIM)
+}
+
+// pub fn dvector_<T>(range: RangeInclusive<T>) -> impl Strategy<Value = DVector<T>>
+// where
+//     RangeInclusive<T>: Strategy<Value = T>,
+//     T: Scalar + PartialEq + Copy,
+//     DefaultAllocator: Allocator<T, Dynamic>,
+// {
+//     vector(range, PROPTEST_MATRIX_DIM)
+// }
+
+macro_rules! define_strategies(
+    ($($strategy_: ident $strategy: ident<$nrows: ident, $ncols: ident>),*) => {$(
+        #[allow(dead_code)]
+        pub fn $strategy() -> impl Strategy<Value = MatrixMN<f64, $nrows, $ncols>> {
+            matrix(PROPTEST_F64, $nrows, $ncols)
+        }
+
+        #[allow(dead_code)]
+        pub fn $strategy_<ScalarStrategy>(scalar_strategy: ScalarStrategy) -> impl Strategy<Value = MatrixMN<ScalarStrategy::Value, $nrows, $ncols>>
+            where
+                ScalarStrategy: Strategy + Clone + 'static,
+                ScalarStrategy::Value: Scalar,
+                DefaultAllocator: Allocator<ScalarStrategy::Value, $nrows, $ncols> {
+            matrix(scalar_strategy, $nrows, $ncols)
+        }
+    )*}
+);
+
+define_strategies!(
+    matrix1_ matrix1<U1, U1>,
+    matrix2_ matrix2<U2, U2>,
+    matrix3_ matrix3<U3, U3>,
+    matrix4_ matrix4<U4, U4>,
+    matrix5_ matrix5<U5, U5>,
+    matrix6_ matrix6<U6, U6>,
+
+    matrix5x2_ matrix5x2<U5, U2>,
+    matrix2x5_ matrix2x5<U2, U5>,
+    matrix5x3_ matrix5x3<U5, U3>,
+    matrix3x5_ matrix3x5<U3, U5>,
+    matrix5x4_ matrix5x4<U5, U4>,
+    matrix4x5_ matrix4x5<U4, U5>,
+
+    vector1_ vector1<U1, U1>,
+    vector2_ vector2<U2, U1>,
+    vector3_ vector3<U3, U1>,
+    vector4_ vector4<U4, U1>,
+    vector5_ vector5<U5, U1>,
+    vector6_ vector6<U6, U1>
+);
 
 /// Generate a proptest that tests that all matrices generated with the
 /// provided rows and columns conform to the constraints defined by the

tests/sparse/cs_conversion.rs

@@ -43,7 +43,6 @@ fn cs_matrix_from_triplet() {
     );
 
     let cs_mat = CsMatrix::from_triplet(4, 5, &irows, &icols, &vals);
-    println!("Mat from triplet: {:?}", cs_mat);
     assert!(cs_mat.is_sorted());
     assert_eq!(cs_mat, cs_expected);
 
@@ -62,7 +61,6 @@ fn cs_matrix_from_triplet() {
     }
 
     let cs_mat = CsMatrix::from_triplet(4, 5, &irows, &icols, &vals);
-    println!("Mat from triplet: {:?}", cs_mat);
     assert!(cs_mat.is_sorted());
     assert_eq!(cs_mat, cs_expected);
 
@@ -80,7 +78,6 @@ fn cs_matrix_from_triplet() {
     vals.append(&mut va);
 
     let cs_mat = CsMatrix::from_triplet(4, 5, &irows, &icols, &vals);
-    println!("Mat from triplet: {:?}", cs_mat);
     assert!(cs_mat.is_sorted());
     assert_eq!(cs_mat, cs_expected * 2.0);
 

tests/sparse/cs_matrix_market.rs

@@ -41,7 +41,6 @@ fn cs_matrix_market() {
 "#;
 
     let cs_mat = io::cs_matrix_from_matrix_market_str(file_str).unwrap();
-    println!("CS mat: {:?}", cs_mat);
     let mat: DMatrix<_> = cs_mat.into();
     let expected = DMatrix::from_row_slice(5, 5, &[
         1.0, 0.0,   0.0,   6.0,    0.0,