diff --git a/rand_distr/src/binomial.rs b/rand_distr/src/binomial.rs
index 0e6bf9a1e18..9c2b110ae5a 100644
--- a/rand_distr/src/binomial.rs
+++ b/rand_distr/src/binomial.rs
@@ -326,4 +326,22 @@ mod test {
     fn test_binomial_invalid_lambda_neg() {
         Binomial::new(20, -10.0).unwrap();
     }
+    
+    #[test]
+    fn value_stability() {
+        fn test_samples(n: u64, p: f64, expected: &[u64]) {
+            let distr = Binomial::new(n, p).unwrap();
+            let mut rng = crate::test::rng(353);
+            let mut buf = [0; 4];
+            for x in &mut buf {
+                *x = rng.sample(&distr);
+            }
+            assert_eq!(buf, expected);
+        }
+        
+        // We have multiple code paths: np < 10, p > 0.5
+        test_samples(2, 0.7, &[1, 1, 2, 1]);
+        test_samples(20, 0.3, &[7, 7, 5, 7]);
+        test_samples(2000, 0.6, &[1194, 1208, 1192, 1210]);
+    }
 }
diff --git a/rand_distr/src/cauchy.rs b/rand_distr/src/cauchy.rs
index 6b0e7c6cf33..8932a8e7787 100644
--- a/rand_distr/src/cauchy.rs
+++ b/rand_distr/src/cauchy.rs
@@ -72,8 +72,7 @@ where Standard: Distribution<N>
 
 #[cfg(test)]
 mod test {
-    use crate::Distribution;
-    use super::Cauchy;
+    use super::*;
 
     fn median(mut numbers: &mut [f64]) -> f64 {
         sort(&mut numbers);
@@ -117,4 +116,22 @@ mod test {
     fn test_cauchy_invalid_scale_neg() {
         Cauchy::new(0.0, -10.0).unwrap();
     }
+    
+    #[test]
+    fn value_stability() {
+        fn test_samples<N: Float + core::fmt::Debug>(m: N, s: N, expected: &[N])
+        where Standard: Distribution<N> {
+            let distr = Cauchy::new(m, s).unwrap();
+            let mut rng = crate::test::rng(353);
+            let mut buf = [m; 4];
+            for x in &mut buf {
+                *x = rng.sample(&distr);
+            }
+            assert_eq!(buf, expected);
+        }
+        
+        test_samples(100f64, 10.0, &[77.93369152808678, 90.1606912098641,
+                125.31516221323625, 86.10217834773925]);
+        test_samples(20f32, 10.0, &[27.175842, -2.0663052, 11.013268, 10.160688]);
+    }
 }
diff --git a/rand_distr/src/dirichlet.rs b/rand_distr/src/dirichlet.rs
index 71cf73c9033..cba063ae755 100644
--- a/rand_distr/src/dirichlet.rs
+++ b/rand_distr/src/dirichlet.rs
@@ -107,8 +107,7 @@ where StandardNormal: Distribution<N>, Exp1: Distribution<N>, Open01: Distributi
 
 #[cfg(test)]
 mod test {
-    use super::Dirichlet;
-    use crate::Distribution;
+    use super::*;
 
     #[test]
     fn test_dirichlet() {
@@ -151,4 +150,14 @@ mod test {
     fn test_dirichlet_invalid_alpha() {
         Dirichlet::new_with_size(0.0f64, 2).unwrap();
     }
+    
+    #[test]
+    fn value_stability() {
+        let mut rng = crate::test::rng(223);
+        assert_eq!(rng.sample(Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap()),
+                vec![0.12941567177708177, 0.4702121891675036, 0.4003721390554146]);
+        assert_eq!(rng.sample(Dirichlet::new_with_size(8.0, 5).unwrap()),
+                vec![0.17684200044809556, 0.29915953935953055,
+                    0.1832858056608014, 0.1425623503573967, 0.19815030417417595]);
+    }
 }
diff --git a/rand_distr/src/exponential.rs b/rand_distr/src/exponential.rs
index 83224890b05..8695252a53f 100644
--- a/rand_distr/src/exponential.rs
+++ b/rand_distr/src/exponential.rs
@@ -121,8 +121,7 @@ where Exp1: Distribution<N>
 
 #[cfg(test)]
 mod test {
-    use crate::Distribution;
-    use super::Exp;
+    use super::*;
 
     #[test]
     fn test_exp() {
@@ -142,4 +141,28 @@ mod test {
     fn test_exp_invalid_lambda_neg() {
         Exp::new(-10.0).unwrap();
     }
+    
+    #[test]
+    fn value_stability() {
+        fn test_samples<N: Float + core::fmt::Debug, D: Distribution<N>>
+        (distr: D, zero: N, expected: &[N])
+        {
+            let mut rng = crate::test::rng(223);
+            let mut buf = [zero; 4];
+            for x in &mut buf {
+                *x = rng.sample(&distr);
+            }
+            assert_eq!(buf, expected);
+        }
+        
+        test_samples(Exp1, 0f32, &[1.079617, 1.8325565, 0.04601166, 0.34471703]);
+        test_samples(Exp1, 0f64, &[1.0796170642388276, 1.8325565304274,
+                0.04601166186842716, 0.3447170217100157]);
+        
+        test_samples(Exp::new(2.0).unwrap(), 0f32,
+                &[0.5398085, 0.91627824, 0.02300583, 0.17235851]);
+        test_samples(Exp::new(1.0).unwrap(), 0f64, &[
+                1.0796170642388276, 1.8325565304274,
+                0.04601166186842716, 0.3447170217100157]);
+    }
 }
diff --git a/rand_distr/src/gamma.rs b/rand_distr/src/gamma.rs
index 4018361648e..bbf861bd537 100644
--- a/rand_distr/src/gamma.rs
+++ b/rand_distr/src/gamma.rs
@@ -417,8 +417,7 @@ where StandardNormal: Distribution<N>, Exp1: Distribution<N>, Open01: Distributi
 
 #[cfg(test)]
 mod test {
-    use crate::Distribution;
-    use super::{Beta, ChiSquared, StudentT, FisherF};
+    use super::*;
 
     #[test]
     fn test_chi_squared_one() {
@@ -482,4 +481,60 @@ mod test {
     fn test_beta_invalid_dof() {
         Beta::new(0., 0.).unwrap();
     }
+    
+    #[test]
+    fn value_stability() {
+        fn test_samples<N: Float + core::fmt::Debug, D: Distribution<N>>
+        (distr: D, zero: N, expected: &[N])
+        {
+            let mut rng = crate::test::rng(223);
+            let mut buf = [zero; 4];
+            for x in &mut buf {
+                *x = rng.sample(&distr);
+            }
+            assert_eq!(buf, expected);
+        }
+        
+        // Gamma has 3 cases: shape == 1, shape < 1, shape > 1
+        test_samples(Gamma::new(1.0, 5.0).unwrap(), 0f32,
+                &[5.398085, 9.162783, 0.2300583, 1.7235851]);
+        test_samples(Gamma::new(0.8, 5.0).unwrap(), 0f32,
+                &[0.5051203, 0.9048302, 3.095812, 1.8566116]);
+        test_samples(Gamma::new(1.1, 5.0).unwrap(), 0f64, &[
+                7.783878094584059, 1.4939528171618057,
+                8.638017638857592, 3.0949337228829004]);
+        
+        // ChiSquared has 2 cases: k == 1, k != 1
+        test_samples(ChiSquared::new(1.0).unwrap(), 0f64, &[
+                0.4893526200348249, 1.635249736808788,
+                0.5013580219361969, 0.1457735613733489]);
+        test_samples(ChiSquared::new(0.1).unwrap(), 0f64, &[
+                0.014824404726978617, 0.021602123937134326,
+                0.0000003431429746851693, 0.00000002291755769542258]);
+        test_samples(ChiSquared::new(10.0).unwrap(), 0f32, 
+                &[12.693656, 6.812016, 11.082001, 12.436167]);
+        
+        // FisherF has same special cases as ChiSquared on each param
+        test_samples(FisherF::new(1.0, 13.5).unwrap(), 0f32,
+                &[0.32283646, 0.048049655, 0.0788893, 1.817178]);
+        test_samples(FisherF::new(1.0, 1.0).unwrap(), 0f32,
+                &[0.29925257, 3.4392934, 9.567652, 0.020074]);
+        test_samples(FisherF::new(0.7, 13.5).unwrap(), 0f64, &[
+                3.3196593155045124, 0.3409169916262829,
+                0.03377989856426519, 0.00004041672861036937]);
+        
+        // StudentT has same special cases as ChiSquared
+        test_samples(StudentT::new(1.0).unwrap(), 0f32,
+                &[0.54703987, -1.8545331, 3.093162, -0.14168274]);
+        test_samples(StudentT::new(1.1).unwrap(), 0f64, &[
+                0.7729195887949754, 1.2606210611616204,
+                -1.7553606501113175, -2.377641221169782]);
+        
+        // Beta has same special cases as Gamma on each param
+        test_samples(Beta::new(1.0, 0.8).unwrap(), 0f32,
+                &[0.6444564, 0.357635, 0.4110078, 0.7347192]);
+        test_samples(Beta::new(0.7, 1.2).unwrap(), 0f64, &[
+                0.6433129944095513, 0.5373371199711573,
+                0.10313293199269491, 0.002472280249144378]);
+    }
 }
diff --git a/rand_distr/src/normal.rs b/rand_distr/src/normal.rs
index 882754fae44..0413ceaf541 100644
--- a/rand_distr/src/normal.rs
+++ b/rand_distr/src/normal.rs
@@ -185,8 +185,7 @@ where StandardNormal: Distribution<N>
 
 #[cfg(test)]
 mod tests {
-    use crate::Distribution;
-    use super::{Normal, LogNormal};
+    use super::*;
 
     #[test]
     fn test_normal() {
@@ -216,4 +215,36 @@ mod tests {
     fn test_log_normal_invalid_sd() {
         LogNormal::new(10.0, -1.0).unwrap();
     }
+    
+    #[test]
+    fn value_stability() {
+        fn test_samples<N: Float + core::fmt::Debug, D: Distribution<N>>
+        (distr: D, zero: N, expected: &[N])
+        {
+            let mut rng = crate::test::rng(213);
+            let mut buf = [zero; 4];
+            for x in &mut buf {
+                *x = rng.sample(&distr);
+            }
+            assert_eq!(buf, expected);
+        }
+        
+        test_samples(StandardNormal, 0f32,
+                &[-0.11844189, 0.781378, 0.06563994, -1.1932899]);
+        test_samples(StandardNormal, 0f64, &[
+                -0.11844188827977231, 0.7813779637772346,
+                0.06563993969580051, -1.1932899004186373]);
+        
+        test_samples(Normal::new(0.0, 1.0).unwrap(), 0f32,
+                &[-0.11844189, 0.781378, 0.06563994, -1.1932899]);
+        test_samples(Normal::new(2.0, 0.5).unwrap(), 0f64, &[
+                1.940779055860114, 2.3906889818886174,
+                2.0328199698479, 1.4033550497906813]);
+        
+        test_samples(LogNormal::new(0.0, 1.0).unwrap(), 0f32,
+                &[0.88830346, 2.1844804, 1.0678421, 0.30322206]);
+        test_samples(LogNormal::new(2.0, 0.5).unwrap(), 0f64, &[
+                6.964174338639032, 10.921015733601452,
+                7.6355881556915906, 4.068828213584092]);
+    }
 }
diff --git a/rand_distr/src/pareto.rs b/rand_distr/src/pareto.rs
index 33ea382d383..ccbce9c2566 100644
--- a/rand_distr/src/pareto.rs
+++ b/rand_distr/src/pareto.rs
@@ -66,8 +66,7 @@ where OpenClosed01: Distribution<N>
 
 #[cfg(test)]
 mod tests {
-    use crate::Distribution;
-    use super::Pareto;
+    use super::*;
 
     #[test]
     #[should_panic]
@@ -86,4 +85,24 @@ mod tests {
             assert!(r >= scale);
         }
     }
+    
+    #[test]
+    fn value_stability() {
+        fn test_samples<N: Float + core::fmt::Debug, D: Distribution<N>>
+        (distr: D, zero: N, expected: &[N])
+        {
+            let mut rng = crate::test::rng(213);
+            let mut buf = [zero; 4];
+            for x in &mut buf {
+                *x = rng.sample(&distr);
+            }
+            assert_eq!(buf, expected);
+        }
+        
+        test_samples(Pareto::new(1.0, 1.0).unwrap(), 0f32,
+                &[1.0423688, 2.1235929, 4.132709, 1.4679428]);
+        test_samples(Pareto::new(2.0, 0.5).unwrap(), 0f64, &[
+                9.019295276219136, 4.3097126018270595,
+                6.837815045397157, 105.8826669383772]);
+    }
 }
diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs
index 4f4a0b7a3d9..d720ab4f413 100644
--- a/rand_distr/src/poisson.rs
+++ b/rand_distr/src/poisson.rs
@@ -134,8 +134,7 @@ where Standard: Distribution<N>
 
 #[cfg(test)]
 mod test {
-    use crate::Distribution;
-    use super::Poisson;
+    use super::*;
 
     #[test]
     fn test_poisson_10() {
@@ -230,4 +229,23 @@ mod test {
     fn test_poisson_invalid_lambda_neg() {
         Poisson::new(-10.0).unwrap();
     }
+    
+    #[test]
+    fn value_stability() {
+        fn test_samples<N: Float + core::fmt::Debug, D: Distribution<N>>
+        (distr: D, zero: N, expected: &[N])
+        {
+            let mut rng = crate::test::rng(223);
+            let mut buf = [zero; 4];
+            for x in &mut buf {
+                *x = rng.sample(&distr);
+            }
+            assert_eq!(buf, expected);
+        }
+        
+        // Special cases: < 12, >= 12
+        test_samples(Poisson::new(7.0).unwrap(), 0f32, &[5.0, 11.0, 6.0, 5.0]);
+        test_samples(Poisson::new(7.0).unwrap(), 0f64, &[9.0, 5.0, 7.0, 6.0]);
+        test_samples(Poisson::new(27.0).unwrap(), 0f32, &[28.0, 32.0, 36.0, 36.0]);
+    }
 }
diff --git a/rand_distr/src/weibull.rs b/rand_distr/src/weibull.rs
index ddde3804d71..388503284fa 100644
--- a/rand_distr/src/weibull.rs
+++ b/rand_distr/src/weibull.rs
@@ -63,8 +63,7 @@ where OpenClosed01: Distribution<N>
 
 #[cfg(test)]
 mod tests {
-    use crate::Distribution;
-    use super::Weibull;
+    use super::*;
 
     #[test]
     #[should_panic]
@@ -83,4 +82,24 @@ mod tests {
             assert!(r >= 0.);
         }
     }
+    
+    #[test]
+    fn value_stability() {
+        fn test_samples<N: Float + core::fmt::Debug, D: Distribution<N>>
+        (distr: D, zero: N, expected: &[N])
+        {
+            let mut rng = crate::test::rng(213);
+            let mut buf = [zero; 4];
+            for x in &mut buf {
+                *x = rng.sample(&distr);
+            }
+            assert_eq!(buf, expected);
+        }
+        
+        test_samples(Weibull::new(1.0, 1.0).unwrap(), 0f32,
+                &[0.041495778, 0.7531094, 1.4189332, 0.38386202]);
+        test_samples(Weibull::new(2.0, 0.5).unwrap(), 0f64, &[
+                1.1343478702739669, 0.29470010050655226,
+                0.7556151370284702, 7.877212340241561]);
+    }
 }