diff --git a/Cargo.toml b/Cargo.toml
index c61c1e71d..7d74b4f38 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name    = "nalgebra"
-version = "0.18.1"
+version = "0.18.2"
 authors = [ "Sébastien Crozet <developer@crozet.re>" ]
 
 description = "Linear algebra library with transformations and statically-sized or dynamically-sized matrices."
@@ -20,32 +20,35 @@ name = "nalgebra"
 path = "src/lib.rs"
 
 [features]
-default         = [ "std" ]
-std             = [ "matrixmultiply", "rand/std", "alga/std" ]
-stdweb          = [ "rand/stdweb" ]
-arbitrary       = [ "quickcheck" ]
+default         = [ "std", "rand_with_std" ]
+std             = [ "matrixmultiply", "alga/std", "alloc" ]
+# Since Cargo doesn't support multiplicitave features, use this to enable rand+std:
+rand_with_std   = [ "rand/std", "std", "rand_distr" ]
+stdweb          = [ "rand", "rand/stdweb" ]
+arbitrary       = [ "quickcheck", "rand" ]
 serde-serialize = [ "serde", "serde_derive", "num-complex/serde" ]
 abomonation-serialize = [ "abomonation" ]
 sparse = [ ]
-debug = [ "approx/num-complex", "rand/std" ]
+debug = [ "rand_with_std" ]
 alloc = [ ]
 io = [ "pest", "pest_derive" ]
+serde_derive = []   # dummy feature to avoid breakage
 
 [dependencies]
 typenum        = "1.10"
 generic-array  = "0.12"
-rand           = { version = "0.6", default-features = false }
+rand           = { version = "0.7", optional = true, default-features = false }
+rand_distr     = { version = "0.2", optional = true }
 num-traits     = { version = "0.2", default-features = false }
 num-complex    = { version = "0.2", default-features = false }
 num-rational   = { version = "0.2", default-features = false }
-approx         = { version = "0.3", default-features = false }
+approx         = { version = "0.3", default-features = false, features = ["num-complex"] }
 alga           = { version = "0.9", default-features = false }
 matrixmultiply = { version = "0.2", optional = true }
-serde          = { version = "1.0", optional = true }
-serde_derive   = { version = "1.0", optional = true }
+serde          = { version = "1.0", optional = true, features = ["derive"] }
 abomonation    = { version = "0.7", optional = true }
 mint           = { version = "0.5", optional = true }
-quickcheck     = { version = "0.8", optional = true }
+quickcheck     = { version = "0.9", optional = true }
 pest           = { version = "2.0", optional = true }
 pest_derive    = { version = "2.0", optional = true }
 
@@ -54,7 +57,9 @@ pest_derive    = { version = "2.0", optional = true }
 
 [dev-dependencies]
 serde_json = "1.0"
-rand_xorshift = "0.1"
+rand = { version = "0.7", features = ["small_rng"] }
+# IsaacRng is used by benches; keep for reproducibility?
+rand_isaac = "0.2"
 ### 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
diff --git a/benches/core/matrix.rs b/benches/core/matrix.rs
index 12fb836d3..ea46e57b8 100644
--- a/benches/core/matrix.rs
+++ b/benches/core/matrix.rs
@@ -1,5 +1,4 @@
 use na::{DMatrix, DVector, Matrix2, Matrix3, Matrix4, MatrixN, Vector2, Vector3, Vector4, U10};
-use rand::{IsaacRng, Rng};
 use std::ops::{Add, Div, Mul, Sub};
 
 #[path = "../common/macros.rs"]
diff --git a/benches/core/vector.rs b/benches/core/vector.rs
index 7d3237e8a..7a4cd8ea5 100644
--- a/benches/core/vector.rs
+++ b/benches/core/vector.rs
@@ -1,5 +1,6 @@
 use na::{DVector, Vector2, Vector3, Vector4, VectorN};
-use rand::{IsaacRng, Rng};
+use rand::{Rng, SeedableRng};
+use rand_isaac::IsaacRng;
 use std::ops::{Add, Div, Mul, Sub};
 use typenum::U10000;
 
@@ -48,7 +49,6 @@ bench_binop_ref!(vec10000_dot_f64, VectorN<f64, U10000>, VectorN<f64, U10000>, d
 bench_binop_ref!(vec10000_dot_f32, VectorN<f32, U10000>, VectorN<f32, U10000>, dot);
 
 fn vec10000_axpy_f64(bh: &mut criterion::Criterion) {
-    use rand::SeedableRng;
     let mut rng = IsaacRng::seed_from_u64(0);
     let mut a = DVector::new_random(10000);
     let b = DVector::new_random(10000);
@@ -58,7 +58,6 @@ fn vec10000_axpy_f64(bh: &mut criterion::Criterion) {
 }
 
 fn vec10000_axpy_beta_f64(bh: &mut criterion::Criterion) {
-    use rand::SeedableRng;
     let mut rng = IsaacRng::seed_from_u64(0);
     let mut a = DVector::new_random(10000);
     let b = DVector::new_random(10000);
@@ -69,7 +68,6 @@ fn vec10000_axpy_beta_f64(bh: &mut criterion::Criterion) {
 }
 
 fn vec10000_axpy_f64_slice(bh: &mut criterion::Criterion) {
-    use rand::SeedableRng;
     let mut rng = IsaacRng::seed_from_u64(0);
     let mut a = DVector::new_random(10000);
     let b = DVector::new_random(10000);
@@ -84,7 +82,6 @@ fn vec10000_axpy_f64_slice(bh: &mut criterion::Criterion) {
 }
 
 fn vec10000_axpy_f64_static(bh: &mut criterion::Criterion) {
-    use rand::SeedableRng;
     let mut rng = IsaacRng::seed_from_u64(0);
     let mut a = VectorN::<f64, U10000>::new_random();
     let b = VectorN::<f64, U10000>::new_random();
@@ -95,7 +92,6 @@ fn vec10000_axpy_f64_static(bh: &mut criterion::Criterion) {
 }
 
 fn vec10000_axpy_f32(bh: &mut criterion::Criterion) {
-    use rand::SeedableRng;
     let mut rng = IsaacRng::seed_from_u64(0);
     let mut a = DVector::new_random(10000);
     let b = DVector::new_random(10000);
@@ -105,7 +101,6 @@ fn vec10000_axpy_f32(bh: &mut criterion::Criterion) {
 }
 
 fn vec10000_axpy_beta_f32(bh: &mut criterion::Criterion) {
-    use rand::SeedableRng;
     let mut rng = IsaacRng::seed_from_u64(0);
     let mut a = DVector::new_random(10000);
     let b = DVector::new_random(10000);
diff --git a/benches/geometry/quaternion.rs b/benches/geometry/quaternion.rs
index dd079aacd..ff9baee11 100644
--- a/benches/geometry/quaternion.rs
+++ b/benches/geometry/quaternion.rs
@@ -1,5 +1,4 @@
 use na::{Quaternion, UnitQuaternion, Vector3};
-use rand::{IsaacRng, Rng};
 use std::ops::{Add, Div, Mul, Sub};
 
 #[path = "../common/macros.rs"]
diff --git a/benches/lib.rs b/benches/lib.rs
index d5333542d..6ca7e8301 100644
--- a/benches/lib.rs
+++ b/benches/lib.rs
@@ -2,23 +2,19 @@
 #![allow(unused_macros)]
 
 extern crate nalgebra as na;
-extern crate rand;
-extern crate test;
-extern crate typenum;
 
 #[macro_use]
 extern crate criterion;
 
 use na::DMatrix;
-use rand::{IsaacRng, Rng};
+use rand::{Rng, SeedableRng};
 
 pub mod core;
 pub mod geometry;
 pub mod linalg;
 
 fn reproductible_dmatrix(nrows: usize, ncols: usize) -> DMatrix<f64> {
-    use rand::SeedableRng;
-    let mut rng = IsaacRng::seed_from_u64(0);
+    let mut rng = rand_isaac::IsaacRng::seed_from_u64(0);
     DMatrix::<f64>::from_fn(nrows, ncols, |_, _| rng.gen())
 }
 
@@ -36,4 +32,4 @@ criterion_main!(
     linalg::solve,
     linalg::svd,
     linalg::symmetric_eigen,
-);
\ No newline at end of file
+);
diff --git a/ci/build.sh b/ci/build.sh
index 550c9a69a..7b36d27c4 100755
--- a/ci/build.sh
+++ b/ci/build.sh
@@ -4,7 +4,9 @@ set -ev
 
 if [ -z "$NO_STD" ]; then
     if [ -z "$LAPACK" ]; then
+        cargo build --verbose -p nalgebra --no-default-features --lib;
         cargo build --verbose -p nalgebra;
+        cargo build --verbose -p nalgebra --features "rand";
         cargo build --verbose -p nalgebra --features "arbitrary";
         cargo build --verbose -p nalgebra --features "mint";
         cargo build --verbose -p nalgebra --features "alloc";
@@ -25,4 +27,4 @@ EOF
     rustup component add rust-src
     cargo install xargo
     xargo build --verbose --no-default-features --target=x86_64-unknown-linux-gnu --features "${CARGO_FEATURES}";
-fi
\ No newline at end of file
+fi
diff --git a/ci/test.sh b/ci/test.sh
index a3a27fb27..ba56b7d87 100755
--- a/ci/test.sh
+++ b/ci/test.sh
@@ -4,11 +4,12 @@ set -ev
 
 if [ -z "$NO_STD" ]; then
     if [ -z "$LAPACK" ]; then
+        cargo test --verbose --no-default-features --features "std" --tests;
         cargo test --verbose;
-        cargo test --verbose "arbitrary";
+        cargo test --verbose --features "arbitrary";
         cargo test --verbose --all-features;
         cd nalgebra-glm; cargo test --verbose;
     else
         cd nalgebra-lapack; cargo test --verbose;
     fi
-fi
\ No newline at end of file
+fi
diff --git a/examples/dimensional_genericity.rs b/examples/dimensional_genericity.rs
index 411a06664..1d96ca744 100644
--- a/examples/dimensional_genericity.rs
+++ b/examples/dimensional_genericity.rs
@@ -1,4 +1,3 @@
-extern crate alga;
 extern crate nalgebra as na;
 
 use alga::linear::FiniteDimInnerSpace;
diff --git a/examples/identity.rs b/examples/identity.rs
index c20c5616c..00f0ec6f5 100644
--- a/examples/identity.rs
+++ b/examples/identity.rs
@@ -1,4 +1,3 @@
-extern crate alga;
 extern crate nalgebra as na;
 
 use alga::linear::Transformation;
@@ -36,4 +35,4 @@ fn main() {
     // They both return the same result.
     assert!(result1 == Vector3::new(100001.0, 200002.0, 300003.0));
     assert!(result2 == Vector3::new(100001.0, 200002.0, 300003.0));
-}
\ No newline at end of file
+}
diff --git a/examples/scalar_genericity.rs b/examples/scalar_genericity.rs
index 75f6f9d40..5ef89248a 100644
--- a/examples/scalar_genericity.rs
+++ b/examples/scalar_genericity.rs
@@ -1,4 +1,3 @@
-extern crate alga;
 extern crate nalgebra as na;
 
 use alga::general::{RealField, RingCommutative};
diff --git a/examples/transform_matrix4.rs b/examples/transform_matrix4.rs
index a9e6cd6e8..74649d74c 100644
--- a/examples/transform_matrix4.rs
+++ b/examples/transform_matrix4.rs
@@ -1,4 +1,3 @@
-extern crate alga;
 #[macro_use]
 extern crate approx;
 extern crate nalgebra as na;
diff --git a/examples/transform_vector_point3.rs b/examples/transform_vector_point3.rs
index 57a189cce..cf2921a5f 100644
--- a/examples/transform_vector_point3.rs
+++ b/examples/transform_vector_point3.rs
@@ -1,4 +1,3 @@
-extern crate alga;
 extern crate nalgebra as na;
 
 use na::{Matrix4, Point3, Vector3, Vector4};
diff --git a/nalgebra-glm/src/geometric.rs b/nalgebra-glm/src/geometric.rs
index ffdc825ae..7c1e1987e 100644
--- a/nalgebra-glm/src/geometric.rs
+++ b/nalgebra-glm/src/geometric.rs
@@ -4,7 +4,7 @@ use crate::aliases::{TVec, TVec3};
 use crate::traits::{Alloc, Dimension, Number};
 
 /// The cross product of two vectors.
-pub fn cross<N: Number, D: Dimension>(x: &TVec3<N>, y: &TVec3<N>) -> TVec3<N> {
+pub fn cross<N: Number>(x: &TVec3<N>, y: &TVec3<N>) -> TVec3<N> {
     x.cross(y)
 }
 
diff --git a/nalgebra-glm/src/lib.rs b/nalgebra-glm/src/lib.rs
index c8f266a45..fe3a32a02 100644
--- a/nalgebra-glm/src/lib.rs
+++ b/nalgebra-glm/src/lib.rs
@@ -116,7 +116,6 @@
 extern crate num_traits as num;
 #[macro_use]
 extern crate approx;
-extern crate alga;
 extern crate nalgebra as na;
 
 pub use crate::aliases::*;
diff --git a/nalgebra-lapack/Cargo.toml b/nalgebra-lapack/Cargo.toml
index 9d3927129..986e8a22c 100644
--- a/nalgebra-lapack/Cargo.toml
+++ b/nalgebra-lapack/Cargo.toml
@@ -13,7 +13,8 @@ license       = "BSD-3-Clause"
 edition       = "2018"
 
 [features]
-serde-serialize = [ "serde", "serde_derive" ]
+serde-serialize = [ "serde" ]   # kept to avoid breakage; prefer to use serde feature directly
+serde_derive = []   # dummy feature to avoid breakage
 
 # For BLAS/LAPACK
 default    = ["openblas"]
@@ -27,14 +28,13 @@ nalgebra      = { version = "0.18", path = ".." }
 num-traits    = "0.2"
 num-complex   = { version = "0.2", default-features = false }
 alga          = { version = "0.9", default-features = false }
-serde         = { version = "1.0", optional = true }
-serde_derive  = { version = "1.0", optional = true }
+serde         = { version = "1.0", optional = true, features = ["derive"] }
 lapack        = { version = "0.16", default-features = false }
 lapack-src    = { version = "0.2", default-features = false }
 # clippy = "*"
 
 [dev-dependencies]
 nalgebra   = { version = "0.18", path = "..", features = [ "arbitrary" ] }
-quickcheck = "0.8"
+quickcheck = "0.9"
 approx     = "0.3"
-rand       = "0.6"
+rand       = "0.7"
diff --git a/nalgebra-lapack/benches/lib.rs b/nalgebra-lapack/benches/lib.rs
index 141e65d98..3b23d2597 100644
--- a/nalgebra-lapack/benches/lib.rs
+++ b/nalgebra-lapack/benches/lib.rs
@@ -2,7 +2,5 @@
 
 extern crate nalgebra as na;
 extern crate nalgebra_lapack as nl;
-extern crate rand;
-extern crate test;
 
 mod linalg;
diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs
index d58aba0bd..39965a59c 100644
--- a/nalgebra-lapack/src/lib.rs
+++ b/nalgebra-lapack/src/lib.rs
@@ -73,11 +73,7 @@
     html_root_url = "https://nalgebra.org/rustdoc"
 )]
 
-extern crate alga;
-extern crate lapack;
-extern crate lapack_src;
 extern crate nalgebra as na;
-extern crate num_complex;
 extern crate num_traits as num;
 
 mod lapack_check;
diff --git a/src/base/alias.rs b/src/base/alias.rs
index a8925cf32..d28af20bc 100644
--- a/src/base/alias.rs
+++ b/src/base/alias.rs
@@ -1,7 +1,7 @@
-#[cfg(any(feature = "alloc", feature = "std"))]
+#[cfg(feature = "alloc")]
 use crate::base::dimension::Dynamic;
 use crate::base::dimension::{U1, U2, U3, U4, U5, U6};
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use crate::base::vec_storage::VecStorage;
 use crate::base::storage::Owned;
 use crate::base::Matrix;
@@ -24,7 +24,7 @@ pub type MatrixMN<N, R, C> = Matrix<N, R, C, Owned<N, R, C>>;
 pub type MatrixN<N, D> = MatrixMN<N, D, D>;
 
 /// A dynamically sized column-major matrix.
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 pub type DMatrix<N> = MatrixN<N, Dynamic>;
 
 /// A stack-allocated, column-major, 1x1 square matrix.
@@ -118,7 +118,7 @@ pub type Matrix6x5<N> = MatrixMN<N, U6, U5>;
  *
  */
 /// A dynamically sized column vector.
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 pub type DVector<N> = Matrix<N, Dynamic, U1, VecStorage<N, Dynamic, U1>>;
 
 /// A statically sized D-dimensional column vector.
@@ -145,7 +145,7 @@ pub type Vector6<N> = VectorN<N, U6>;
  *
  */
 /// A dynamically sized row vector.
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 pub type RowDVector<N> = Matrix<N, U1, Dynamic, VecStorage<N, U1, Dynamic>>;
 
 /// A statically sized D-dimensional row vector.
diff --git a/src/base/blas.rs b/src/base/blas.rs
index beb83c4ba..ded53f7f3 100644
--- a/src/base/blas.rs
+++ b/src/base/blas.rs
@@ -21,8 +21,6 @@ impl<N: ComplexField, D: Dim, S: Storage<N, D>> Vector<N, D, S> {
     /// # Examples:
     ///
     /// ```
-    /// # extern crate num_complex;
-    /// # extern crate nalgebra;
     /// # use num_complex::Complex;
     /// # use nalgebra::Vector3;
     /// let vec = Vector3::new(Complex::new(11.0, 3.0), Complex::new(-15.0, 0.0), Complex::new(13.0, 5.0));
@@ -565,6 +563,14 @@ where
         );
 
         if ncols2 == 0 {
+            // NOTE: we can't just always multiply by beta
+            // because we documented the guaranty that `self` is
+            // never read if `beta` is zero.
+            if beta.is_zero() {
+                self.fill(N::zero());
+            } else {
+                *self *= beta;
+            }
             return;
         }
 
@@ -991,92 +997,109 @@ where N: Scalar + Zero + ClosedAdd + ClosedMul
 
         #[cfg(feature = "std")]
         {
-            // matrixmultiply can be used only if the std feature is available.
-            let nrows1 = self.nrows();
-            let (nrows2, ncols2) = a.shape();
-            let (nrows3, ncols3) = b.shape();
-
-            assert_eq!(
-                ncols2, nrows3,
-                "gemm: dimensions mismatch for multiplication."
-            );
-            assert_eq!(
-                (nrows1, ncols1),
-                (nrows2, ncols3),
-                "gemm: dimensions mismatch for addition."
-            );
-
             // We assume large matrices will be Dynamic but small matrices static.
             // We could use matrixmultiply for large statically-sized matrices but the performance
             // threshold to activate it would be different from SMALL_DIM because our code optimizes
             // better for statically-sized matrices.
-            let is_dynamic = R1::is::<Dynamic>()
+            if R1::is::<Dynamic>()
                 || C1::is::<Dynamic>()
                 || R2::is::<Dynamic>()
                 || C2::is::<Dynamic>()
                 || R3::is::<Dynamic>()
-                || C3::is::<Dynamic>();
-            // Threshold determined empirically.
-            const SMALL_DIM: usize = 5;
-
-            if is_dynamic
-                && nrows1 > SMALL_DIM
-                && ncols1 > SMALL_DIM
-                && nrows2 > SMALL_DIM
-                && ncols2 > SMALL_DIM
-            {
-                if N::is::<f32>() {
-                    let (rsa, csa) = a.strides();
-                    let (rsb, csb) = b.strides();
-                    let (rsc, csc) = self.strides();
-
-                    unsafe {
-                        matrixmultiply::sgemm(
-                            nrows2,
-                            ncols2,
-                            ncols3,
-                            mem::transmute_copy(&alpha),
-                            a.data.ptr() as *const f32,
-                            rsa as isize,
-                            csa as isize,
-                            b.data.ptr() as *const f32,
-                            rsb as isize,
-                            csb as isize,
-                            mem::transmute_copy(&beta),
-                            self.data.ptr_mut() as *mut f32,
-                            rsc as isize,
-                            csc as isize,
-                        );
+                || C3::is::<Dynamic>() {
+                // matrixmultiply can be used only if the std feature is available.
+                let nrows1 = self.nrows();
+                let (nrows2, ncols2) = a.shape();
+                let (nrows3, ncols3) = b.shape();
+
+                // Threshold determined empirically.
+                const SMALL_DIM: usize = 5;
+
+                if nrows1 > SMALL_DIM
+                    && ncols1 > SMALL_DIM
+                    && nrows2 > SMALL_DIM
+                    && ncols2 > SMALL_DIM
+                {
+                    assert_eq!(
+                        ncols2, nrows3,
+                        "gemm: dimensions mismatch for multiplication."
+                    );
+                    assert_eq!(
+                        (nrows1, ncols1),
+                        (nrows2, ncols3),
+                        "gemm: dimensions mismatch for addition."
+                    );
+
+                    // NOTE: this case should never happen because we enter this
+                    // codepath only when ncols2 > SMALL_DIM. Though we keep this
+                    // here just in case if in the future we change the conditions to
+                    // enter this codepath.
+                    if ncols2 == 0 {
+                        // NOTE: we can't just always multiply by beta
+                        // because we documented the guaranty that `self` is
+                        // never read if `beta` is zero.
+                        if beta.is_zero() {
+                            self.fill(N::zero());
+                        } else {
+                            *self *= beta;
+                        }
+                        return;
                     }
-                    return;
-                } else if N::is::<f64>() {
-                    let (rsa, csa) = a.strides();
-                    let (rsb, csb) = b.strides();
-                    let (rsc, csc) = self.strides();
-
-                    unsafe {
-                        matrixmultiply::dgemm(
-                            nrows2,
-                            ncols2,
-                            ncols3,
-                            mem::transmute_copy(&alpha),
-                            a.data.ptr() as *const f64,
-                            rsa as isize,
-                            csa as isize,
-                            b.data.ptr() as *const f64,
-                            rsb as isize,
-                            csb as isize,
-                            mem::transmute_copy(&beta),
-                            self.data.ptr_mut() as *mut f64,
-                            rsc as isize,
-                            csc as isize,
-                        );
+
+                    if N::is::<f32>() {
+                        let (rsa, csa) = a.strides();
+                        let (rsb, csb) = b.strides();
+                        let (rsc, csc) = self.strides();
+
+                        unsafe {
+                            matrixmultiply::sgemm(
+                                nrows2,
+                                ncols2,
+                                ncols3,
+                                mem::transmute_copy(&alpha),
+                                a.data.ptr() as *const f32,
+                                rsa as isize,
+                                csa as isize,
+                                b.data.ptr() as *const f32,
+                                rsb as isize,
+                                csb as isize,
+                                mem::transmute_copy(&beta),
+                                self.data.ptr_mut() as *mut f32,
+                                rsc as isize,
+                                csc as isize,
+                            );
+                        }
+                        return;
+                    } else if N::is::<f64>() {
+                        let (rsa, csa) = a.strides();
+                        let (rsb, csb) = b.strides();
+                        let (rsc, csc) = self.strides();
+
+                        unsafe {
+                            matrixmultiply::dgemm(
+                                nrows2,
+                                ncols2,
+                                ncols3,
+                                mem::transmute_copy(&alpha),
+                                a.data.ptr() as *const f64,
+                                rsa as isize,
+                                csa as isize,
+                                b.data.ptr() as *const f64,
+                                rsb as isize,
+                                csb as isize,
+                                mem::transmute_copy(&beta),
+                                self.data.ptr_mut() as *mut f64,
+                                rsc as isize,
+                                csc as isize,
+                            );
+                        }
+                        return;
                     }
-                    return;
                 }
             }
         }
 
+
         for j1 in 0..ncols1 {
             // FIXME: avoid bound checks.
             self.column_mut(j1).gemv(alpha, a, &b.column(j1), beta);
diff --git a/src/base/construction.rs b/src/base/construction.rs
index 5eace4bb0..7ed3bef0e 100644
--- a/src/base/construction.rs
+++ b/src/base/construction.rs
@@ -4,14 +4,12 @@ use crate::base::storage::Owned;
 use quickcheck::{Arbitrary, Gen};
 
 use num::{Bounded, One, Zero};
-use rand::distributions::{Distribution, Standard};
-use rand::Rng;
-#[cfg(feature = "std")]
-use rand::{self, distributions::StandardNormal};
+#[cfg(feature = "rand")]
+use rand::{Rng, distributions::{Distribution, Standard}};
 use std::iter;
 use typenum::{self, Cmp, Greater};
 
-#[cfg(feature = "std")]
+#[cfg(feature = "rand_distr")]
 use alga::general::RealField;
 use alga::general::{ClosedAdd, ClosedMul};
 
@@ -240,13 +238,15 @@ where DefaultAllocator: Allocator<N, R, C>
 
     /// Creates a matrix filled with random values.
     #[inline]
-    #[cfg(feature = "std")]
+    #[cfg(feature = "rand_with_std")]
     pub fn new_random_generic(nrows: R, ncols: C) -> Self
     where Standard: Distribution<N> {
-        Self::from_fn_generic(nrows, ncols, |_, _| rand::random())
+        let mut rng = rand::thread_rng();
+        Self::from_fn_generic(nrows, ncols, |_, _| rng.gen())
     }
 
     /// Creates a matrix filled with random values from the given distribution.
+    #[cfg(feature = "rand")]
     #[inline]
     pub fn from_distribution_generic<Distr: Distribution<N> + ?Sized, G: Rng + ?Sized>(
         nrows: R,
@@ -547,6 +547,7 @@ macro_rules! impl_constructors(
             }
 
             /// Creates a matrix or vector filled with random values from the given distribution.
+            #[cfg(feature = "rand")]
             #[inline]
             pub fn from_distribution<Distr: Distribution<N> + ?Sized, G: Rng + ?Sized>(
                 $($args: usize,)*
@@ -557,6 +558,7 @@ macro_rules! impl_constructors(
             }
         }
 
+        #[cfg(feature = "rand")]
         impl<N: Scalar, $($DimIdent: $DimBound, )*> MatrixMN<N $(, $Dims)*>
             where
             DefaultAllocator: Allocator<N $(, $Dims)*>,
@@ -564,7 +566,7 @@ macro_rules! impl_constructors(
 
             /// Creates a matrix filled with random values.
             #[inline]
-            #[cfg(feature = "std")]
+            #[cfg(feature = "rand_with_std")]
             pub fn new_random($($args: usize),*) -> Self {
                 Self::new_random_generic($($gargs),*)
             }
@@ -760,6 +762,7 @@ where
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: Scalar, R: Dim, C: Dim> Distribution<MatrixMN<N, R, C>> for Standard
 where
     DefaultAllocator: Allocator<N, R, C>,
@@ -774,7 +777,7 @@ where
     }
 }
 
-#[cfg(feature = "arbitrary")]
+#[cfg(all(feature = "arbitrary", feature = "rand"))]
 impl<N, R, C> Arbitrary for MatrixMN<N, R, C>
 where
     R: Dim,
@@ -794,11 +797,12 @@ where
     }
 }
 
-#[cfg(feature = "std")]
+// TODO(specialization): faster impls possible for D≤4 (see rand_distr::{UnitCircle, UnitSphere})
+#[cfg(feature = "rand_distr")]
 impl<N: RealField, D: DimName> Distribution<Unit<VectorN<N, D>>> for Standard
 where
     DefaultAllocator: Allocator<N, D>,
-    StandardNormal: Distribution<N>,
+    rand_distr::StandardNormal: Distribution<N>,
 {
     /// Generate a uniformly distributed random unit vector.
     #[inline]
@@ -806,7 +810,7 @@ where
         Unit::new_normalize(VectorN::from_distribution_generic(
             D::name(),
             U1,
-            &StandardNormal,
+            &rand_distr::StandardNormal,
             rng,
         ))
     }
diff --git a/src/base/conversion.rs b/src/base/conversion.rs
index 4c5bb0170..93a49fbb0 100644
--- a/src/base/conversion.rs
+++ b/src/base/conversion.rs
@@ -14,11 +14,11 @@ use crate::base::constraint::{SameNumberOfColumns, SameNumberOfRows, ShapeConstr
 use crate::base::dimension::{
     Dim, DimName, U1, U10, U11, U12, U13, U14, U15, U16, U2, U3, U4, U5, U6, U7, U8, U9,
 };
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use crate::base::dimension::Dynamic;
 use crate::base::iter::{MatrixIter, MatrixIterMut};
 use crate::base::storage::{ContiguousStorage, ContiguousStorageMut, Storage, StorageMut};
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use crate::base::VecStorage;
 use crate::base::{DefaultAllocator, Matrix, ArrayStorage, MatrixMN, MatrixSlice, MatrixSliceMut, Scalar};
 
@@ -353,7 +353,7 @@ where
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<'a, N, C, RStride, CStride> From<MatrixSlice<'a, N, Dynamic, C, RStride, CStride>>
     for Matrix<N, Dynamic, C, VecStorage<N, Dynamic, C>>
 where
@@ -367,7 +367,7 @@ where
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<'a, N, R, RStride, CStride> From<MatrixSlice<'a, N, R, Dynamic, RStride, CStride>>
     for Matrix<N, R, Dynamic, VecStorage<N, R, Dynamic>>
 where
@@ -397,7 +397,7 @@ where
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<'a, N, C, RStride, CStride> From<MatrixSliceMut<'a, N, Dynamic, C, RStride, CStride>>
     for Matrix<N, Dynamic, C, VecStorage<N, Dynamic, C>>
 where
@@ -411,7 +411,7 @@ where
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<'a, N, R, RStride, CStride> From<MatrixSliceMut<'a, N, R, Dynamic, RStride, CStride>>
     for Matrix<N, R, Dynamic, VecStorage<N, R, Dynamic>>
 where
diff --git a/src/base/coordinates.rs b/src/base/coordinates.rs
index 832723e39..be0929cfc 100644
--- a/src/base/coordinates.rs
+++ b/src/base/coordinates.rs
@@ -6,6 +6,8 @@
 
 use std::mem;
 use std::ops::{Deref, DerefMut};
+#[cfg(feature = "serde-serialize")]
+use serde::{Serialize, Deserialize};
 
 use crate::base::dimension::{U1, U2, U3, U4, U5, U6};
 use crate::base::storage::{ContiguousStorage, ContiguousStorageMut};
diff --git a/src/base/default_allocator.rs b/src/base/default_allocator.rs
index c07c87085..a04651aa3 100644
--- a/src/base/default_allocator.rs
+++ b/src/base/default_allocator.rs
@@ -15,11 +15,11 @@ use generic_array::ArrayLength;
 use typenum::Prod;
 
 use crate::base::allocator::{Allocator, Reallocator};
-#[cfg(any(feature = "alloc", feature = "std"))]
+#[cfg(feature = "alloc")]
 use crate::base::dimension::Dynamic;
 use crate::base::dimension::{Dim, DimName};
 use crate::base::array_storage::ArrayStorage;
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use crate::base::vec_storage::VecStorage;
 use crate::base::storage::{Storage, StorageMut};
 use crate::base::Scalar;
@@ -75,7 +75,7 @@ where
 
 // Dynamic - Static
 // Dynamic - Dynamic
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, C: Dim> Allocator<N, Dynamic, C> for DefaultAllocator {
     type Buffer = VecStorage<N, Dynamic, C>;
 
@@ -106,7 +106,7 @@ impl<N: Scalar, C: Dim> Allocator<N, Dynamic, C> for DefaultAllocator {
 }
 
 // Static - Dynamic
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, R: DimName> Allocator<N, R, Dynamic> for DefaultAllocator {
     type Buffer = VecStorage<N, R, Dynamic>;
 
@@ -172,7 +172,7 @@ where
 }
 
 // Static × Static -> Dynamic × Any
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, RFrom, CFrom, CTo> Reallocator<N, RFrom, CFrom, Dynamic, CTo> for DefaultAllocator
 where
     RFrom: DimName,
@@ -201,7 +201,7 @@ where
 }
 
 // Static × Static -> Static × Dynamic
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, RFrom, CFrom, RTo> Reallocator<N, RFrom, CFrom, RTo, Dynamic> for DefaultAllocator
 where
     RFrom: DimName,
@@ -230,7 +230,7 @@ where
 }
 
 // All conversion from a dynamic buffer to a dynamic buffer.
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, CFrom: Dim, CTo: Dim> Reallocator<N, Dynamic, CFrom, Dynamic, CTo>
     for DefaultAllocator
 {
@@ -246,7 +246,7 @@ impl<N: Scalar, CFrom: Dim, CTo: Dim> Reallocator<N, Dynamic, CFrom, Dynamic, CT
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, CFrom: Dim, RTo: DimName> Reallocator<N, Dynamic, CFrom, RTo, Dynamic>
     for DefaultAllocator
 {
@@ -262,7 +262,7 @@ impl<N: Scalar, CFrom: Dim, RTo: DimName> Reallocator<N, Dynamic, CFrom, RTo, Dy
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, RFrom: DimName, CTo: Dim> Reallocator<N, RFrom, Dynamic, Dynamic, CTo>
     for DefaultAllocator
 {
@@ -278,7 +278,7 @@ impl<N: Scalar, RFrom: DimName, CTo: Dim> Reallocator<N, RFrom, Dynamic, Dynamic
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, RFrom: DimName, RTo: DimName> Reallocator<N, RFrom, Dynamic, RTo, Dynamic>
     for DefaultAllocator
 {
diff --git a/src/base/edition.rs b/src/base/edition.rs
index e473246fd..6db834907 100644
--- a/src/base/edition.rs
+++ b/src/base/edition.rs
@@ -1,20 +1,20 @@
 use num::{One, Zero};
 use std::cmp;
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use std::iter::ExactSizeIterator;
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use std::mem;
 use std::ptr;
 
 use crate::base::allocator::{Allocator, Reallocator};
 use crate::base::constraint::{DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint};
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use crate::base::dimension::Dynamic;
 use crate::base::dimension::{
     Dim, DimAdd, DimDiff, DimMin, DimMinimum, DimName, DimSub, DimSum, U1,
 };
 use crate::base::storage::{Storage, StorageMut};
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use crate::base::DMatrix;
 use crate::base::{DefaultAllocator, Matrix, MatrixMN, RowVector, Scalar, Vector};
 
@@ -40,7 +40,7 @@ impl<N: Scalar + Zero, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     }
 
     /// Creates a new matrix by extracting the given set of rows from `self`.
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn select_rows<'a, I>(&self, irows: I) -> MatrixMN<N, Dynamic, C>
     where
         I: IntoIterator<Item = &'a usize>,
@@ -72,7 +72,7 @@ impl<N: Scalar + Zero, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     }
 
     /// Creates a new matrix by extracting the given set of columns from `self`.
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn select_columns<'a, I>(&self, icols: I) -> MatrixMN<N, R, Dynamic>
     where
         I: IntoIterator<Item = &'a usize>,
@@ -308,7 +308,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     }
 
     /// Removes all columns in `indices`   
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn remove_columns_at(self, indices: &[usize]) -> MatrixMN<N, R, Dynamic>
     where
         C: DimSub<Dynamic, Output = Dynamic>,
@@ -345,7 +345,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     }
 
     /// Removes all rows in `indices`   
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn remove_rows_at(self, indices: &[usize]) -> MatrixMN<N, Dynamic, C>
     where
         R: DimSub<Dynamic, Output = Dynamic>,
@@ -392,7 +392,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
 
     /// Removes `n` consecutive columns from this matrix, starting with the `i`-th (included).
     #[inline]
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn remove_columns(self, i: usize, n: usize) -> MatrixMN<N, R, Dynamic>
     where
         C: DimSub<Dynamic, Output = Dynamic>,
@@ -475,7 +475,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
 
     /// Removes `n` consecutive rows from this matrix, starting with the `i`-th (included).
     #[inline]
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn remove_rows(self, i: usize, n: usize) -> MatrixMN<N, Dynamic, C>
     where
         R: DimSub<Dynamic, Output = Dynamic>,
@@ -553,7 +553,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
 
     /// Inserts `n` columns filled with `val` starting at the `i-th` position.
     #[inline]
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn insert_columns(self, i: usize, n: usize, val: N) -> MatrixMN<N, R, Dynamic>
     where
         C: DimAdd<Dynamic, Output = Dynamic>,
@@ -631,7 +631,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
 
     /// Inserts `n` rows filled with `val` starting at the `i-th` position.
     #[inline]
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn insert_rows(self, i: usize, n: usize, val: N) -> MatrixMN<N, Dynamic, C>
     where
         R: DimAdd<Dynamic, Output = Dynamic>,
@@ -691,7 +691,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     ///
     /// The values are copied such that `self[(i, j)] == result[(i, j)]`. If the result has more
     /// rows and/or columns than `self`, then the extra rows or columns are filled with `val`.
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn resize(self, new_nrows: usize, new_ncols: usize, val: N) -> DMatrix<N>
     where DefaultAllocator: Reallocator<N, R, C, Dynamic, Dynamic> {
         self.resize_generic(Dynamic::new(new_nrows), Dynamic::new(new_ncols), val)
@@ -701,7 +701,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     ///
     /// The values are copied such that `self[(i, j)] == result[(i, j)]`. If the result has more
     /// rows than `self`, then the extra rows are filled with `val`.
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn resize_vertically(self, new_nrows: usize, val: N) -> MatrixMN<N, Dynamic, C>
     where DefaultAllocator: Reallocator<N, R, C, Dynamic, C> {
         let ncols = self.data.shape().1;
@@ -712,7 +712,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     ///
     /// The values are copied such that `self[(i, j)] == result[(i, j)]`. If the result has more
     /// columns than `self`, then the extra columns are filled with `val`.
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn resize_horizontally(self, new_ncols: usize, val: N) -> MatrixMN<N, R, Dynamic>
     where DefaultAllocator: Reallocator<N, R, C, R, Dynamic> {
         let nrows = self.data.shape().0;
@@ -796,7 +796,7 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar> DMatrix<N> {
     /// Resizes this matrix in-place.
     ///
@@ -813,7 +813,7 @@ impl<N: Scalar> DMatrix<N> {
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, C: Dim> MatrixMN<N, Dynamic, C>
 where DefaultAllocator: Allocator<N, Dynamic, C>
 {
@@ -823,7 +823,7 @@ where DefaultAllocator: Allocator<N, Dynamic, C>
     /// rows than `self`, then the extra rows are filled with `val`.
     ///
     /// Defined only for owned matrices with a dynamic number of rows (for example, `DVector`).
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn resize_vertically_mut(&mut self, new_nrows: usize, val: N)
     where DefaultAllocator: Reallocator<N, Dynamic, C, Dynamic, C> {
         let placeholder =
@@ -834,7 +834,7 @@ where DefaultAllocator: Allocator<N, Dynamic, C>
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N: Scalar, R: Dim> MatrixMN<N, R, Dynamic>
 where DefaultAllocator: Allocator<N, R, Dynamic>
 {
@@ -844,7 +844,7 @@ where DefaultAllocator: Allocator<N, R, Dynamic>
     /// columns than `self`, then the extra columns are filled with `val`.
     ///
     /// Defined only for owned matrices with a dynamic number of columns (for example, `DVector`).
-    #[cfg(any(feature = "std", feature = "alloc"))]
+    #[cfg(feature = "alloc")]
     pub fn resize_horizontally_mut(&mut self, new_ncols: usize, val: N)
     where DefaultAllocator: Reallocator<N, R, Dynamic, R, Dynamic> {
         let placeholder =
@@ -935,7 +935,7 @@ unsafe fn extend_rows<N: Scalar>(
 
 /// Extend the number of columns of the `Matrix` with elements from
 /// a given iterator.
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N, R, S> Extend<N> for Matrix<N, R, Dynamic, S>
 where
     N: Scalar,
@@ -983,7 +983,7 @@ where
 
 /// Extend the number of rows of the `Vector` with elements from
 /// a given iterator.
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N, S> Extend<N> for Matrix<N, Dynamic, U1, S>
 where
     N: Scalar,
@@ -1004,7 +1004,7 @@ where
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl<N, R, S, RV, SV> Extend<Vector<N, RV, SV>> for Matrix<N, R, Dynamic, S>
 where
     N: Scalar,
diff --git a/src/base/helper.rs b/src/base/helper.rs
index ef85a477b..04a82bb83 100644
--- a/src/base/helper.rs
+++ b/src/base/helper.rs
@@ -1,3 +1,4 @@
+#![cfg(feature = "rand")]
 #[cfg(feature = "arbitrary")]
 use quickcheck::{Arbitrary, Gen};
 use rand::distributions::{Distribution, Standard};
diff --git a/src/base/matrix.rs b/src/base/matrix.rs
index 7119a8ad9..542414361 100644
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@@ -1359,75 +1359,99 @@ where
     }
 }
 
-impl<N, R: Dim, C: Dim, S> fmt::Display for Matrix<N, R, C, S>
-where
-    N: Scalar + fmt::Display,
-    S: Storage<N, R, C>,
-    DefaultAllocator: Allocator<usize, R, C>,
-{
-    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-        #[cfg(feature = "std")]
-        fn val_width<N: Scalar + fmt::Display>(val: N, f: &mut fmt::Formatter) -> usize {
-            match f.precision() {
-                Some(precision) => format!("{:.1$}", val, precision).chars().count(),
-                None => format!("{}", val).chars().count(),
-            }
-        }
-
-        #[cfg(not(feature = "std"))]
-        fn val_width<N: Scalar + fmt::Display>(_: N, _: &mut fmt::Formatter) -> usize {
-            4
-        }
-
-        let (nrows, ncols) = self.data.shape();
+macro_rules! impl_fmt {
+    ($trait: path, $fmt_str_without_precision: expr, $fmt_str_with_precision: expr) => {
+        impl<N, R: Dim, C: Dim, S> $trait for Matrix<N, R, C, S>
+        where
+            N: Scalar + $trait,
+            S: Storage<N, R, C>,
+            DefaultAllocator: Allocator<usize, R, C>,
+        {
+            fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+                #[cfg(feature = "std")]
+                fn val_width<N: Scalar + $trait>(val: N, f: &mut fmt::Formatter) -> usize {
+                    match f.precision() {
+                        Some(precision) => format!($fmt_str_with_precision, val, precision).chars().count(),
+                        None => format!($fmt_str_without_precision, val).chars().count(),
+                    }
+                }
 
-        if nrows.value() == 0 || ncols.value() == 0 {
-            return write!(f, "[ ]");
-        }
+                #[cfg(not(feature = "std"))]
+                fn val_width<N: Scalar + $trait>(_: N, _: &mut fmt::Formatter) -> usize {
+                    4
+                }
 
-        let mut max_length = 0;
-        let mut lengths: MatrixMN<usize, R, C> = Matrix::zeros_generic(nrows, ncols);
-        let (nrows, ncols) = self.shape();
+                let (nrows, ncols) = self.data.shape();
 
-        for i in 0..nrows {
-            for j in 0..ncols {
-                lengths[(i, j)] = val_width(self[(i, j)], f);
-                max_length = crate::max(max_length, lengths[(i, j)]);
-            }
-        }
+                if nrows.value() == 0 || ncols.value() == 0 {
+                    return write!(f, "[ ]");
+                }
 
-        let max_length_with_space = max_length + 1;
+                let mut max_length = 0;
+                let mut lengths: MatrixMN<usize, R, C> = Matrix::zeros_generic(nrows, ncols);
+                let (nrows, ncols) = self.shape();
 
-        writeln!(f)?;
-        writeln!(
-            f,
-            "  ┌ {:>width$} ┐",
-            "",
-            width = max_length_with_space * ncols - 1
-        )?;
+                for i in 0..nrows {
+                    for j in 0..ncols {
+                        lengths[(i, j)] = val_width(self[(i, j)], f);
+                        max_length = crate::max(max_length, lengths[(i, j)]);
+                    }
+                }
 
-        for i in 0..nrows {
-            write!(f, "  │")?;
-            for j in 0..ncols {
-                let number_length = lengths[(i, j)] + 1;
-                let pad = max_length_with_space - number_length;
-                write!(f, " {:>thepad$}", "", thepad = pad)?;
-                match f.precision() {
-                    Some(precision) => write!(f, "{:.1$}", (*self)[(i, j)], precision)?,
-                    None => write!(f, "{}", (*self)[(i, j)])?,
+                let max_length_with_space = max_length + 1;
+
+                writeln!(f)?;
+                writeln!(
+                    f,
+                    "  ┌ {:>width$} ┐",
+                    "",
+                    width = max_length_with_space * ncols - 1
+                )?;
+
+                for i in 0..nrows {
+                    write!(f, "  │")?;
+                    for j in 0..ncols {
+                        let number_length = lengths[(i, j)] + 1;
+                        let pad = max_length_with_space - number_length;
+                        write!(f, " {:>thepad$}", "", thepad = pad)?;
+                        match f.precision() {
+                            Some(precision) => write!(f, $fmt_str_with_precision, (*self)[(i, j)], precision)?,
+                            None => write!(f, $fmt_str_without_precision, (*self)[(i, j)])?,
+                        }
+                    }
+                    writeln!(f, " │")?;
                 }
+
+                writeln!(
+                    f,
+                    "  └ {:>width$} ┘",
+                    "",
+                    width = max_length_with_space * ncols - 1
+                )?;
+                writeln!(f)
             }
-            writeln!(f, " │")?;
         }
-
-        writeln!(
-            f,
-            "  └ {:>width$} ┘",
-            "",
-            width = max_length_with_space * ncols - 1
-        )?;
-        writeln!(f)
-    }
+    };
+}
+impl_fmt!(fmt::Display, "{}", "{:.1$}");
+impl_fmt!(fmt::LowerExp, "{:e}", "{:.1$e}");
+impl_fmt!(fmt::UpperExp, "{:E}", "{:.1$E}");
+impl_fmt!(fmt::Octal, "{:o}", "{:1$o}");
+impl_fmt!(fmt::LowerHex, "{:x}", "{:1$x}");
+impl_fmt!(fmt::UpperHex, "{:X}", "{:1$X}");
+impl_fmt!(fmt::Binary, "{:b}", "{:.1$b}");
+impl_fmt!(fmt::Pointer, "{:p}", "{:.1$p}");
+
+#[test]
+fn lower_exp() {
+    let test = crate::Matrix2::new(1e6, 2e5, 2e-5, 1.);
+    assert_eq!(format!("{:e}", test), r"
+  ┌           ┐
+  │  1e6  2e5 │
+  │ 2e-5  1e0 │
+  └           ┘
+
+")
 }
 
 impl<N: Scalar + Ring, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
diff --git a/src/base/matrix_alga.rs b/src/base/matrix_alga.rs
index ac6aced78..5419f02d6 100644
--- a/src/base/matrix_alga.rs
+++ b/src/base/matrix_alga.rs
@@ -285,7 +285,7 @@ where DefaultAllocator: Allocator<N, R, C>
                 }
             }
             _ => {
-                #[cfg(any(feature = "std", feature = "alloc"))]
+                #[cfg(feature = "alloc")]
                 {
                     // XXX: use a GenericArray instead.
                     let mut known_basis = Vec::new();
@@ -310,7 +310,7 @@ where DefaultAllocator: Allocator<N, R, C>
                         }
                     }
                 }
-                #[cfg(all(not(feature = "std"), not(feature = "alloc")))]
+                #[cfg(not(feature = "alloc"))]
                 {
                     panic!("Cannot compute the orthogonal subspace basis of a vector with a dimension greater than 3 \
                             if #![no_std] is enabled and the 'alloc' feature is not enabled.")
diff --git a/src/base/mod.rs b/src/base/mod.rs
index 0ec6311c7..139e560ba 100644
--- a/src/base/mod.rs
+++ b/src/base/mod.rs
@@ -23,7 +23,7 @@ mod matrix;
 mod matrix_alga;
 mod array_storage;
 mod matrix_slice;
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 mod vec_storage;
 mod properties;
 mod scalar;
@@ -47,5 +47,5 @@ pub use self::alias::*;
 pub use self::alias_slice::*;
 pub use self::array_storage::*;
 pub use self::matrix_slice::*;
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 pub use self::vec_storage::*;
diff --git a/src/base/ops.rs b/src/base/ops.rs
index f389eae52..ef0014c5c 100644
--- a/src/base/ops.rs
+++ b/src/base/ops.rs
@@ -918,10 +918,11 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     /// # use nalgebra::Vector3;
     /// assert_eq!(Vector3::new(-1.0, 2.0, 3.0).max(), 3.0);
     /// assert_eq!(Vector3::new(-1.0, -2.0, -3.0).max(), -1.0);
+    /// assert_eq!(Vector3::new(5u32, 2, 3).max(), 5);
     /// ```
     #[inline]
     pub fn max(&self) -> N
-        where N: PartialOrd + Signed {
+        where N: PartialOrd + Zero {
         self.xcmp(|e| e, Ordering::Greater)
     }
 
@@ -959,10 +960,11 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
     /// # use nalgebra::Vector3;
     /// assert_eq!(Vector3::new(-1.0, 2.0, 3.0).min(), -1.0);
     /// assert_eq!(Vector3::new(1.0, 2.0, 3.0).min(), 1.0);
+    /// assert_eq!(Vector3::new(5u32, 2, 3).min(), 2);
     /// ```
     #[inline]
     pub fn min(&self) -> N
-        where N: PartialOrd + Signed {
+        where N: PartialOrd + Zero {
         self.xcmp(|e| e, Ordering::Less)
     }
 }
\ No newline at end of file
diff --git a/src/base/vec_storage.rs b/src/base/vec_storage.rs
index 2b4bf7432..a8f7f28c7 100644
--- a/src/base/vec_storage.rs
+++ b/src/base/vec_storage.rs
@@ -4,6 +4,9 @@ use std::io::{Result as IOResult, Write};
 #[cfg(all(feature = "alloc", not(feature = "std")))]
 use alloc::vec::Vec;
 
+#[cfg(feature = "serde-serialize")]
+use serde::{Serialize, Deserialize};
+
 use crate::base::allocator::Allocator;
 use crate::base::default_allocator::DefaultAllocator;
 use crate::base::dimension::{Dim, DimName, Dynamic, U1};
diff --git a/src/geometry/isometry_construction.rs b/src/geometry/isometry_construction.rs
index 979c39556..25562d0a7 100644
--- a/src/geometry/isometry_construction.rs
+++ b/src/geometry/isometry_construction.rs
@@ -4,8 +4,8 @@ use crate::base::storage::Owned;
 use quickcheck::{Arbitrary, Gen};
 
 use num::One;
-use rand::distributions::{Distribution, Standard};
-use rand::Rng;
+#[cfg(feature = "rand")]
+use rand::{Rng, distributions::{Distribution, Standard}};
 
 use alga::general::RealField;
 use alga::linear::Rotation as AlgaRotation;
@@ -75,6 +75,7 @@ where DefaultAllocator: Allocator<N, D>
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: RealField, D: DimName, R> Distribution<Isometry<N, D, R>> for Standard
 where
     R: AlgaRotation<Point<N, D>>,
diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs
index 2c7131182..0976ad97d 100644
--- a/src/geometry/orthographic.rs
+++ b/src/geometry/orthographic.rs
@@ -1,7 +1,7 @@
 #[cfg(feature = "arbitrary")]
 use quickcheck::{Arbitrary, Gen};
-use rand::distributions::{Distribution, Standard};
-use rand::Rng;
+#[cfg(feature = "rand")]
+use rand::{Rng, distributions::{Distribution, Standard}};
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 use std::fmt;
@@ -10,6 +10,7 @@ use std::mem;
 use alga::general::RealField;
 
 use crate::base::dimension::U3;
+#[cfg(feature = "rand")]
 use crate::base::helper;
 use crate::base::storage::Storage;
 use crate::base::{Matrix4, Vector, Vector3};
@@ -678,9 +679,11 @@ impl<N: RealField> Orthographic3<N> {
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: RealField> Distribution<Orthographic3<N>> for Standard
 where Standard: Distribution<N>
 {
+    /// Generate an arbitrary random variate for testing purposes.
     fn sample<R: Rng + ?Sized>(&self, r: &mut R) -> Orthographic3<N> {
         let left = r.gen();
         let right = helper::reject_rand(r, |x: &N| *x > left);
diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs
index 8020c0cf5..cfc173cad 100644
--- a/src/geometry/perspective.rs
+++ b/src/geometry/perspective.rs
@@ -1,7 +1,7 @@
 #[cfg(feature = "arbitrary")]
 use quickcheck::{Arbitrary, Gen};
-use rand::distributions::{Distribution, Standard};
-use rand::Rng;
+#[cfg(feature = "rand")]
+use rand::{Rng, distributions::{Distribution, Standard}};
 
 #[cfg(feature = "serde-serialize")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
@@ -11,6 +11,7 @@ use std::mem;
 use alga::general::RealField;
 
 use crate::base::dimension::U3;
+#[cfg(feature = "rand")]
 use crate::base::helper;
 use crate::base::storage::Storage;
 use crate::base::{Matrix4, Scalar, Vector, Vector3};
@@ -261,9 +262,11 @@ impl<N: RealField> Perspective3<N> {
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: RealField> Distribution<Perspective3<N>> for Standard
 where Standard: Distribution<N>
 {
+    /// Generate an arbitrary random variate for testing purposes.
     fn sample<'a, R: Rng + ?Sized>(&self, r: &'a mut R) -> Perspective3<N> {
         let znear = r.gen();
         let zfar = helper::reject_rand(r, |&x: &N| !(x - znear).is_zero());
diff --git a/src/geometry/point_construction.rs b/src/geometry/point_construction.rs
index 2fac11d42..33b97078c 100644
--- a/src/geometry/point_construction.rs
+++ b/src/geometry/point_construction.rs
@@ -2,8 +2,8 @@
 use quickcheck::{Arbitrary, Gen};
 
 use num::{Bounded, One, Zero};
-use rand::distributions::{Distribution, Standard};
-use rand::Rng;
+#[cfg(feature = "rand")]
+use rand::{Rng, distributions::{Distribution, Standard}};
 
 use alga::general::ClosedDiv;
 use crate::base::allocator::Allocator;
@@ -126,11 +126,13 @@ where DefaultAllocator: Allocator<N, D>
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: Scalar, D: DimName> Distribution<Point<N, D>> for Standard
 where
     DefaultAllocator: Allocator<N, D>,
     Standard: Distribution<N>,
 {
+    /// Generates a `Point` where each coordinate is an independent variate from `[0, 1)`.
     #[inline]
     fn sample<'a, G: Rng + ?Sized>(&self, rng: &mut G) -> Point<N, D> {
         Point::from(rng.gen::<VectorN<N, D>>())
diff --git a/src/geometry/quaternion_construction.rs b/src/geometry/quaternion_construction.rs
index 893fdb6ae..387cc5d07 100644
--- a/src/geometry/quaternion_construction.rs
+++ b/src/geometry/quaternion_construction.rs
@@ -6,8 +6,8 @@ use crate::base::storage::Owned;
 use quickcheck::{Arbitrary, Gen};
 
 use num::{One, Zero};
-use rand::distributions::{Distribution, OpenClosed01, Standard};
-use rand::Rng;
+#[cfg(feature = "rand")]
+use rand::{Rng, distributions::{Distribution, OpenClosed01, Standard}};
 
 use alga::general::RealField;
 
@@ -124,6 +124,7 @@ impl<N: RealField> Zero for Quaternion<N> {
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: RealField> Distribution<Quaternion<N>> for Standard
 where Standard: Distribution<N>
 {
@@ -685,6 +686,7 @@ impl<N: RealField> One for UnitQuaternion<N> {
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: RealField> Distribution<UnitQuaternion<N>> for Standard
 where OpenClosed01: Distribution<N>
 {
@@ -722,15 +724,14 @@ where
     }
 }
 
-#[cfg(test)]
+#[cfg(all(test, feature="rand"))]
 mod tests {
-    extern crate rand_xorshift;
     use super::*;
-    use rand::SeedableRng;
+    use rand::{Rng, SeedableRng, rngs::SmallRng};
 
     #[test]
     fn random_unit_quats_are_unit() {
-        let mut rng = rand_xorshift::XorShiftRng::from_seed([0xAB; 16]);
+        let mut rng = SmallRng::seed_from_u64(2);
         for _ in 0..1000 {
             let x = rng.gen::<UnitQuaternion<f32>>();
             assert!(relative_eq!(x.into_inner().norm(), 1.0))
diff --git a/src/geometry/rotation_specialization.rs b/src/geometry/rotation_specialization.rs
index 4baec27cb..e6a4530fb 100644
--- a/src/geometry/rotation_specialization.rs
+++ b/src/geometry/rotation_specialization.rs
@@ -5,7 +5,9 @@ use quickcheck::{Arbitrary, Gen};
 
 use alga::general::RealField;
 use num::Zero;
-use rand::distributions::{Distribution, OpenClosed01, Standard};
+#[cfg(feature = "rand")]
+use rand::distributions::{Distribution, OpenClosed01, Standard, uniform::SampleUniform};
+#[cfg(feature = "rand")]
 use rand::Rng;
 use std::ops::Neg;
 
@@ -230,13 +232,14 @@ impl<N: RealField> Rotation2<N> {
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: RealField> Distribution<Rotation2<N>> for Standard
-where OpenClosed01: Distribution<N>
+where N: SampleUniform
 {
     /// Generate a uniformly distributed random rotation.
     #[inline]
     fn sample<'a, R: Rng + ?Sized>(&self, rng: &'a mut R) -> Rotation2<N> {
-        Rotation2::new(rng.sample(OpenClosed01) * N::two_pi())
+        Rotation2::new(rng.gen_range(N::zero(), N::two_pi()))
     }
 }
 
@@ -817,8 +820,9 @@ impl<N: RealField> Rotation3<N> {
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: RealField> Distribution<Rotation3<N>> for Standard
-where OpenClosed01: Distribution<N>
+where OpenClosed01: Distribution<N>, N: SampleUniform
 {
     /// Generate a uniformly distributed random rotation.
     #[inline]
@@ -828,7 +832,7 @@ where OpenClosed01: Distribution<N>
         // In D. Kirk, editor, Graphics Gems III, pages 117-120. Academic, New York, 1992.
 
         // Compute a random rotation around Z
-        let theta = N::two_pi() * rng.sample(OpenClosed01);
+        let theta = rng.gen_range(N::zero(), N::two_pi());
         let (ts, tc) = theta.sin_cos();
         let a = MatrixN::<N, U3>::new(
             tc,
@@ -843,7 +847,7 @@ where OpenClosed01: Distribution<N>
         );
 
         // Compute a random rotation *of* Z
-        let phi = N::two_pi() * rng.sample(OpenClosed01);
+        let phi = rng.gen_range(N::zero(), N::two_pi());
         let z = rng.sample(OpenClosed01);
         let (ps, pc) = phi.sin_cos();
         let sqrt_z = z.sqrt();
diff --git a/src/geometry/similarity_construction.rs b/src/geometry/similarity_construction.rs
index a3722ba95..959bd67e5 100644
--- a/src/geometry/similarity_construction.rs
+++ b/src/geometry/similarity_construction.rs
@@ -4,8 +4,8 @@ use crate::base::storage::Owned;
 use quickcheck::{Arbitrary, Gen};
 
 use num::One;
-use rand::distributions::{Distribution, Standard};
-use rand::Rng;
+#[cfg(feature = "rand")]
+use rand::{Rng, distributions::{Distribution, Standard}};
 
 use alga::general::RealField;
 use alga::linear::Rotation as AlgaRotation;
@@ -57,12 +57,14 @@ where
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: RealField, D: DimName, R> Distribution<Similarity<N, D, R>> for Standard
 where
     R: AlgaRotation<Point<N, D>>,
     DefaultAllocator: Allocator<N, D>,
     Standard: Distribution<N> + Distribution<R>,
 {
+    /// Generate an arbitrary random variate for testing purposes.
     #[inline]
     fn sample<'a, G: Rng + ?Sized>(&self, rng: &mut G) -> Similarity<N, D, R> {
         let mut s = rng.gen();
diff --git a/src/geometry/translation_construction.rs b/src/geometry/translation_construction.rs
index 339bdd2a6..e2285befa 100644
--- a/src/geometry/translation_construction.rs
+++ b/src/geometry/translation_construction.rs
@@ -4,8 +4,8 @@ use crate::base::storage::Owned;
 use quickcheck::{Arbitrary, Gen};
 
 use num::{One, Zero};
-use rand::distributions::{Distribution, Standard};
-use rand::Rng;
+#[cfg(feature = "rand")]
+use rand::{Rng, distributions::{Distribution, Standard}};
 
 use alga::general::ClosedAdd;
 
@@ -47,11 +47,13 @@ where DefaultAllocator: Allocator<N, D>
     }
 }
 
+#[cfg(feature = "rand")]
 impl<N: Scalar, D: DimName> Distribution<Translation<N, D>> for Standard
 where
     DefaultAllocator: Allocator<N, D>,
     Standard: Distribution<N>,
 {
+    /// Generate an arbitrary random variate for testing purposes.
     #[inline]
     fn sample<'a, G: Rng + ?Sized>(&self, rng: &'a mut G) -> Translation<N, D> {
         Translation::from(rng.gen::<VectorN<N, D>>())
diff --git a/src/geometry/unit_complex.rs b/src/geometry/unit_complex.rs
index 7ba7f3740..434355ded 100755
--- a/src/geometry/unit_complex.rs
+++ b/src/geometry/unit_complex.rs
@@ -84,7 +84,6 @@ impl<N: RealField> UnitComplex<N> {
     ///
     /// # Example
     /// ```
-    /// # extern crate num_complex;
     /// # use num_complex::Complex;
     /// # use nalgebra::UnitComplex;
     /// let angle = 1.78f32;
@@ -373,4 +372,4 @@ impl<N: RealField> UlpsEq for UnitComplex<N> {
         self.re.ulps_eq(&other.re, epsilon, max_ulps)
             && self.im.ulps_eq(&other.im, epsilon, max_ulps)
     }
-}
\ No newline at end of file
+}
diff --git a/src/geometry/unit_complex_construction.rs b/src/geometry/unit_complex_construction.rs
index 046142c8a..2934816c2 100644
--- a/src/geometry/unit_complex_construction.rs
+++ b/src/geometry/unit_complex_construction.rs
@@ -3,7 +3,9 @@ use quickcheck::{Arbitrary, Gen};
 
 use num::One;
 use num_complex::Complex;
-use rand::distributions::{Distribution, OpenClosed01, Standard};
+#[cfg(feature = "rand_distr")]
+use rand::distributions::{Distribution, Standard};
+#[cfg(feature = "rand_distr")]
 use rand::Rng;
 
 use alga::general::RealField;
@@ -275,13 +277,15 @@ impl<N: RealField> One for UnitComplex<N> {
     }
 }
 
+#[cfg(feature = "rand_distr")]
 impl<N: RealField> Distribution<UnitComplex<N>> for Standard
-where OpenClosed01: Distribution<N>
+where rand_distr::UnitCircle: Distribution<[N; 2]>
 {
     /// Generate a uniformly distributed random `UnitComplex`.
     #[inline]
     fn sample<'a, R: Rng + ?Sized>(&self, rng: &mut R) -> UnitComplex<N> {
-        UnitComplex::from_angle(rng.sample(OpenClosed01) * N::two_pi())
+        let x = rand_distr::UnitCircle.sample(rng);
+        UnitComplex::new_unchecked(Complex::new(x[0], x[1]))
     }
 }
 
diff --git a/src/lib.rs b/src/lib.rs
index a92e2250b..139b53d2b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -82,41 +82,15 @@ an optimized set of tools for computer graphics and physics. Those features incl
 #![deny(unused_qualifications)]
 #![deny(unused_results)]
 #![deny(missing_docs)]
-#![warn(incoherent_fundamental_impls)]
 #![doc(
     html_favicon_url = "https://nalgebra.org/img/favicon.ico",
     html_root_url = "https://nalgebra.org/rustdoc"
 )]
 #![cfg_attr(not(feature = "std"), no_std)]
-#![cfg_attr(all(feature = "alloc", not(feature = "std")), feature(alloc))]
-
-#[cfg(feature = "arbitrary")]
-extern crate quickcheck;
-
-#[cfg(feature = "serde")]
-extern crate serde;
-#[cfg(feature = "serde")]
-#[macro_use]
-extern crate serde_derive;
-
-#[cfg(feature = "abomonation-serialize")]
-extern crate abomonation;
-
-#[cfg(feature = "mint")]
-extern crate mint;
 
 #[macro_use]
 extern crate approx;
-extern crate generic_array;
-#[cfg(feature = "std")]
-extern crate matrixmultiply;
-extern crate num_complex;
 extern crate num_traits as num;
-extern crate num_rational;
-extern crate rand;
-extern crate typenum;
-
-extern crate alga;
 
 #[cfg(all(feature = "alloc", not(feature = "std")))]
 extern crate alloc;
@@ -124,8 +98,6 @@ extern crate alloc;
 #[cfg(not(feature = "std"))]
 extern crate core as std;
 
-#[cfg(feature = "io")]
-extern crate pest;
 #[macro_use]
 #[cfg(feature = "io")]
 extern crate pest_derive;
diff --git a/src/linalg/permutation_sequence.rs b/src/linalg/permutation_sequence.rs
index ce493905e..3202e0b28 100644
--- a/src/linalg/permutation_sequence.rs
+++ b/src/linalg/permutation_sequence.rs
@@ -6,7 +6,7 @@ use num::One;
 
 use crate::allocator::Allocator;
 use crate::base::{DefaultAllocator, Matrix, Scalar, VectorN};
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 use crate::dimension::Dynamic;
 use crate::dimension::{Dim, DimName, U1};
 use crate::storage::StorageMut;
@@ -51,7 +51,7 @@ where DefaultAllocator: Allocator<(usize, usize), D>
     }
 }
 
-#[cfg(any(feature = "std", feature = "alloc"))]
+#[cfg(feature = "alloc")]
 impl PermutationSequence<Dynamic>
 where DefaultAllocator: Allocator<(usize, usize), Dynamic>
 {
diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs
index fc493dca2..729b2b7c8 100644
--- a/src/linalg/symmetric_eigen.rs
+++ b/src/linalg/symmetric_eigen.rs
@@ -354,6 +354,7 @@ mod test {
         }
     }
 
+    #[cfg(feature = "rand_with_std")]
     #[test]
     fn wilkinson_shift_random() {
         for _ in 0..1000 {
diff --git a/tests/core/empty.rs b/tests/core/empty.rs
new file mode 100644
index 000000000..3e17ee1bb
--- /dev/null
+++ b/tests/core/empty.rs
@@ -0,0 +1,60 @@
+use na::{DMatrix, DVector};
+
+#[test]
+fn empty_matrix_mul_vector() {
+    // Issue #644
+    let m = DMatrix::<f32>::zeros(8, 0);
+    let v = DVector::<f32>::zeros(0);
+    assert_eq!(m * v, DVector::zeros(8));
+}
+
+#[test]
+fn empty_matrix_mul_matrix() {
+    let m1 = DMatrix::<f32>::zeros(3, 0);
+    let m2 = DMatrix::<f32>::zeros(0, 4);
+    assert_eq!(m1 * m2, DMatrix::zeros(3, 4));
+
+    // Still works with larger matrices.
+    let m1 = DMatrix::<f32>::zeros(13, 0);
+    let m2 = DMatrix::<f32>::zeros(0, 14);
+    assert_eq!(m1 * m2, DMatrix::zeros(13, 14));
+}
+
+#[test]
+fn empty_matrix_tr_mul_vector() {
+    let m = DMatrix::<f32>::zeros(0, 5);
+    let v = DVector::<f32>::zeros(0);
+    assert_eq!(m.tr_mul(&v), DVector::zeros(5));
+}
+
+#[test]
+fn empty_matrix_tr_mul_matrix() {
+    let m1 = DMatrix::<f32>::zeros(0, 3);
+    let m2 = DMatrix::<f32>::zeros(0, 4);
+    assert_eq!(m1.tr_mul(&m2), DMatrix::zeros(3, 4));
+}
+
+#[test]
+fn empty_matrix_gemm() {
+    let mut res = DMatrix::repeat(3, 4, 1.0);
+    let m1 = DMatrix::<f32>::zeros(3, 0);
+    let m2 = DMatrix::<f32>::zeros(0, 4);
+    res.gemm(1.0, &m1, &m2, 0.5);
+    assert_eq!(res, DMatrix::repeat(3, 4, 0.5));
+
+    // Still works with lager matrices.
+    let mut res = DMatrix::repeat(13, 14, 1.0);
+    let m1 = DMatrix::<f32>::zeros(13, 0);
+    let m2 = DMatrix::<f32>::zeros(0, 14);
+    res.gemm(1.0, &m1, &m2, 0.5);
+    assert_eq!(res, DMatrix::repeat(13, 14, 0.5));
+}
+
+#[test]
+fn empty_matrix_gemm_tr() {
+    let mut res = DMatrix::repeat(3, 4, 1.0);
+    let m1 = DMatrix::<f32>::zeros(0, 3);
+    let m2 = DMatrix::<f32>::zeros(0, 4);
+    res.gemm_tr(1.0, &m1, &m2, 0.5);
+    assert_eq!(res, DMatrix::repeat(3, 4, 0.5));
+}
\ No newline at end of file
diff --git a/tests/core/helper.rs b/tests/core/helper.rs
index 625a4a464..bbece5695 100644
--- a/tests/core/helper.rs
+++ b/tests/core/helper.rs
@@ -33,6 +33,7 @@ impl<N: RealField> Distribution<RandComplex<N>> for Standard
 
 // This is a wrapper similar to RandComplex, but for non-complex.
 // This exists only to make generic tests easier to write.
+// Generates variates in the range [0, 1). Do we want this? E.g. we could use standard normal samples instead.
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
 pub struct RandScalar<N>(pub N);
 
diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs
index e4fb4d0c2..712802a3d 100644
--- a/tests/core/matrix.rs
+++ b/tests/core/matrix.rs
@@ -466,6 +466,7 @@ fn zip_map() {
     assert_eq!(computed, expected);
 }
 
+#[cfg(feature = "rand_with_std")]
 #[test]
 #[should_panic]
 fn trace_panic() {
diff --git a/tests/core/mod.rs b/tests/core/mod.rs
index c53493bdc..f8c7e4773 100644
--- a/tests/core/mod.rs
+++ b/tests/core/mod.rs
@@ -8,6 +8,7 @@ mod matrix_slice;
 #[cfg(feature = "mint")]
 mod mint;
 mod serde;
+mod empty;
 
 
 #[cfg(feature = "arbitrary")]
diff --git a/tests/lib.rs b/tests/lib.rs
index f6634d0fc..4aa957447 100644
--- a/tests/lib.rs
+++ b/tests/lib.rs
@@ -1,18 +1,10 @@
-#[cfg(feature = "abomonation-serialize")]
-extern crate abomonation;
-extern crate alga;
 #[macro_use]
 extern crate approx;
-#[cfg(feature = "mint")]
-extern crate mint;
 extern crate nalgebra as na;
 extern crate num_traits as num;
 #[cfg(feature = "arbitrary")]
 #[macro_use]
 extern crate quickcheck;
-extern crate rand;
-extern crate serde_json;
-extern crate num_complex;
 
 mod core;
 mod geometry;
diff --git a/tests/linalg/inverse.rs b/tests/linalg/inverse.rs
index 61f433246..f0be4dd72 100644
--- a/tests/linalg/inverse.rs
+++ b/tests/linalg/inverse.rs
@@ -130,7 +130,7 @@ fn matrix5_try_inverse_scaled_identity() {
                                             0.0, 1.0e+20,     0.0,     0.0,     0.0,
                                             0.0,     0.0, 1.0e+20,     0.0,     0.0,
                                             0.0,     0.0,     0.0, 1.0e+20,     0.0,
-                                            0.0,     0.0,     0.0,     0.0, 1.0e+20);;
+                                            0.0,     0.0,     0.0,     0.0, 1.0e+20);
     let a_inv = a.try_inverse().expect("Matrix should be invertible");
 
     assert_relative_eq!(a_inv, expected_inverse);