diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs
index ba071dfa..4682948c 100644
--- a/src/distribution/beta.rs
+++ b/src/distribution/beta.rs
@@ -3,7 +3,6 @@ use crate::function::{beta, gamma};
 use crate::is_zero;
 use crate::statistics::*;
 use crate::{Result, StatsError};
-use core::f64::INFINITY as INF;
 use rand::Rng;
 
 /// Implements the [Beta](https://en.wikipedia.org/wiki/Beta_distribution)
@@ -326,13 +325,13 @@ impl Continuous<f64, f64> for Beta {
             0.0
         } else if self.shape_a.is_infinite() {
             if ulps_eq!(x, 1.0) {
-                INF
+                f64::INFINITY
             } else {
                 0.0
             }
         } else if self.shape_b.is_infinite() {
             if is_zero(x) {
-                INF
+                f64::INFINITY
             } else {
                 0.0
             }
@@ -361,18 +360,18 @@ impl Continuous<f64, f64> for Beta {
     /// where `α` is shapeA, `β` is shapeB, and `Γ` is the gamma function
     fn ln_pdf(&self, x: f64) -> f64 {
         if !(0.0..=1.0).contains(&x) {
-            -INF
+            f64::NEG_INFINITY
         } else if self.shape_a.is_infinite() {
             if ulps_eq!(x, 1.0) {
-                INF
+                f64::INFINITY
             } else {
-                -INF
+                f64::NEG_INFINITY
             }
         } else if self.shape_b.is_infinite() {
             if is_zero(x) {
-                INF
+                f64::INFINITY
             } else {
-                -INF
+                f64::NEG_INFINITY
             }
         } else if ulps_eq!(self.shape_a, 1.0) && ulps_eq!(self.shape_b, 1.0) {
             0.0
@@ -383,14 +382,14 @@ impl Continuous<f64, f64> for Beta {
             let bb = if ulps_eq!(self.shape_a, 1.0) && is_zero(x) {
                 0.0
             } else if is_zero(x) {
-                -INF
+                f64::NEG_INFINITY
             } else {
                 (self.shape_a - 1.0) * x.ln()
             };
             let cc = if ulps_eq!(self.shape_b, 1.0) && ulps_eq!(x, 1.0) {
                 0.0
             } else if ulps_eq!(x, 1.0) {
-                -INF
+                f64::NEG_INFINITY
             } else {
                 (self.shape_b - 1.0) * (1.0 - x).ln()
             };
@@ -412,7 +411,7 @@ mod tests {
 
     #[test]
     fn test_create() {
-        let valid = [(1.0, 1.0), (9.0, 1.0), (5.0, 100.0), (1.0, INF), (INF, 1.0)];
+        let valid = [(1.0, 1.0), (9.0, 1.0), (5.0, 100.0), (1.0, f64::INFINITY), (f64::INFINITY, 1.0)];
         for &arg in valid.iter() {
             try_create(arg);
         }
@@ -424,15 +423,15 @@ mod tests {
             (0.0, 0.0),
             (0.0, 0.1),
             (1.0, 0.0),
-            (0.0, INF),
-            (INF, 0.0),
+            (0.0, f64::INFINITY),
+            (f64::INFINITY, 0.0),
             (f64::NAN, 1.0),
             (1.0, f64::NAN),
             (f64::NAN, f64::NAN),
             (1.0, -1.0),
             (-1.0, 1.0),
             (-1.0, -1.0),
-            (INF, INF),
+            (f64::INFINITY, f64::INFINITY),
         ];
         for &arg in invalid.iter() {
             bad_create_case(arg);
@@ -446,8 +445,8 @@ mod tests {
             ((1.0, 1.0), 0.5),
             ((9.0, 1.0), 0.9),
             ((5.0, 100.0), 0.047619047619047619047616),
-            ((1.0, INF), 0.0),
-            ((INF, 1.0), 1.0),
+            ((1.0, f64::INFINITY), 0.0),
+            ((f64::INFINITY, 1.0), 1.0),
         ];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
@@ -461,8 +460,8 @@ mod tests {
             ((1.0, 1.0), 1.0 / 12.0),
             ((9.0, 1.0), 9.0 / 1100.0),
             ((5.0, 100.0), 500.0 / 1168650.0),
-            ((1.0, INF), 0.0),
-            ((INF, 1.0), 0.0),
+            ((1.0, f64::INFINITY), 0.0),
+            ((f64::INFINITY, 1.0), 0.0),
         ];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
@@ -481,8 +480,8 @@ mod tests {
         }
         test_case_special((1.0, 1.0), 0.0, 1e-14, f);
         let entropy = |x: Beta| x.entropy();
-        test_none((1.0, INF), entropy);
-        test_none((INF, 1.0), entropy);
+        test_none((1.0, f64::INFINITY), entropy);
+        test_none((f64::INFINITY, 1.0), entropy);
     }
 
     #[test]
@@ -491,16 +490,16 @@ mod tests {
         test_case((1.0, 1.0), 0.0, skewness);
         test_case((9.0, 1.0), -1.4740554623801777107177478829, skewness);
         test_case((5.0, 100.0), 0.817594109275534303545831591, skewness);
-        test_case((1.0, INF), 2.0, skewness);
-        test_case((INF, 1.0), -2.0, skewness);
+        test_case((1.0, f64::INFINITY), 2.0, skewness);
+        test_case((f64::INFINITY, 1.0), -2.0, skewness);
     }
 
     #[test]
     fn test_mode() {
         let mode = |x: Beta| x.mode().unwrap();
         test_case((5.0, 100.0), 0.038834951456310676243255386, mode);
-        test_case((92.0, INF), 0.0, mode);
-        test_case((INF, 2.0), 1.0, mode);
+        test_case((92.0, f64::INFINITY), 0.0, mode);
+        test_case((f64::INFINITY, 2.0), 1.0, mode);
     }
 
     #[test]
@@ -539,12 +538,12 @@ mod tests {
             ((5.0, 100.0), 0.5, 4.534102298350337661e-23),
             ((5.0, 100.0), 1.0, 0.0),
             ((5.0, 100.0), 1.0, 0.0),
-            ((1.0, INF), 0.0, INF),
-            ((1.0, INF), 0.5, 0.0),
-            ((1.0, INF), 1.0, 0.0),
-            ((INF, 1.0), 0.0, 0.0),
-            ((INF, 1.0), 0.5, 0.0),
-            ((INF, 1.0), 1.0, INF),
+            ((1.0, f64::INFINITY), 0.0, f64::INFINITY),
+            ((1.0, f64::INFINITY), 0.5, 0.0),
+            ((1.0, f64::INFINITY), 1.0, 0.0),
+            ((f64::INFINITY, 1.0), 0.0, 0.0),
+            ((f64::INFINITY, 1.0), 0.5, 0.0),
+            ((f64::INFINITY, 1.0), 1.0, f64::INFINITY),
         ];
         for &(arg, x, expect) in test.iter() {
             test_case(arg, expect, f(x));
@@ -570,18 +569,18 @@ mod tests {
             ((1.0, 1.0), 0.0, 0.0),
             ((1.0, 1.0), 0.5, 0.0),
             ((1.0, 1.0), 1.0, 0.0),
-            ((9.0, 1.0), 0.0, -INF),
+            ((9.0, 1.0), 0.0, f64::NEG_INFINITY),
             ((9.0, 1.0), 0.5, -3.347952867143343092547366497),
             ((9.0, 1.0), 1.0, 2.1972245773362193827904904738),
-            ((5.0, 100.0), 0.0, -INF),
+            ((5.0, 100.0), 0.0, f64::NEG_INFINITY),
             ((5.0, 100.0), 0.5, -51.447830024537682154565870),
-            ((5.0, 100.0), 1.0, -INF),
-            ((1.0, INF), 0.0, INF),
-            ((1.0, INF), 0.5, -INF),
-            ((1.0, INF), 1.0, -INF),
-            ((INF, 1.0), 0.0, -INF),
-            ((INF, 1.0), 0.5, -INF),
-            ((INF, 1.0), 1.0, INF),
+            ((5.0, 100.0), 1.0, f64::NEG_INFINITY),
+            ((1.0, f64::INFINITY), 0.0, f64::INFINITY),
+            ((1.0, f64::INFINITY), 0.5, f64::NEG_INFINITY),
+            ((1.0, f64::INFINITY), 1.0, f64::NEG_INFINITY),
+            ((f64::INFINITY, 1.0), 0.0, f64::NEG_INFINITY),
+            ((f64::INFINITY, 1.0), 0.5, f64::NEG_INFINITY),
+            ((f64::INFINITY, 1.0), 1.0, f64::INFINITY),
         ];
         for &(arg, x, expect) in test.iter() {
             test_case(arg, expect, f(x));
@@ -591,13 +590,13 @@ mod tests {
     #[test]
     fn test_ln_pdf_input_lt_0() {
         let ln_pdf = |arg: f64| move |x: Beta| x.ln_pdf(arg);
-        test_case((1.0, 1.0), -INF, ln_pdf(-1.0));
+        test_case((1.0, 1.0), f64::NEG_INFINITY, ln_pdf(-1.0));
     }
 
     #[test]
     fn test_ln_pdf_input_gt_1() {
         let ln_pdf = |arg: f64| move |x: Beta| x.ln_pdf(arg);
-        test_case((1.0, 1.0), -INF, ln_pdf(2.0));
+        test_case((1.0, 1.0), f64::NEG_INFINITY, ln_pdf(2.0));
     }
 
     #[test]
@@ -613,12 +612,12 @@ mod tests {
             ((5.0, 100.0), 0.0, 0.0),
             ((5.0, 100.0), 0.5, 1.0),
             ((5.0, 100.0), 1.0, 1.0),
-            ((1.0, INF), 0.0, 1.0),
-            ((1.0, INF), 0.5, 1.0),
-            ((1.0, INF), 1.0, 1.0),
-            ((INF, 1.0), 0.0, 0.0),
-            ((INF, 1.0), 0.5, 0.0),
-            ((INF, 1.0), 1.0, 1.0),
+            ((1.0, f64::INFINITY), 0.0, 1.0),
+            ((1.0, f64::INFINITY), 0.5, 1.0),
+            ((1.0, f64::INFINITY), 1.0, 1.0),
+            ((f64::INFINITY, 1.0), 0.0, 0.0),
+            ((f64::INFINITY, 1.0), 0.5, 0.0),
+            ((f64::INFINITY, 1.0), 1.0, 1.0),
         ];
         for &(arg, x, expect) in test.iter() {
             test_case(arg, expect, cdf(x));
@@ -638,12 +637,12 @@ mod tests {
             ((5.0, 100.0), 0.0, 1.0),
             ((5.0, 100.0), 0.5, 0.0),
             ((5.0, 100.0), 1.0, 0.0),
-            ((1.0, INF), 0.0, 0.0),
-            ((1.0, INF), 0.5, 0.0),
-            ((1.0, INF), 1.0, 0.0),
-            ((INF, 1.0), 0.0, 1.0),
-            ((INF, 1.0), 0.5, 1.0),
-            ((INF, 1.0), 1.0, 0.0),
+            ((1.0, f64::INFINITY), 0.0, 0.0),
+            ((1.0, f64::INFINITY), 0.5, 0.0),
+            ((1.0, f64::INFINITY), 1.0, 0.0),
+            ((f64::INFINITY, 1.0), 0.0, 1.0),
+            ((f64::INFINITY, 1.0), 0.5, 1.0),
+            ((f64::INFINITY, 1.0), 1.0, 0.0),
         ];
         for &(arg, x, expect) in test.iter() {
             test_case(arg, expect, sf(x));
diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs
index 7c5c9622..f3bfc76d 100644
--- a/src/distribution/binomial.rs
+++ b/src/distribution/binomial.rs
@@ -49,7 +49,7 @@ impl Binomial {
     /// assert!(result.is_err());
     /// ```
     pub fn new(p: f64, n: u64) -> Result<Binomial> {
-        if p.is_nan() || p < 0.0 || p > 1.0 {
+        if p.is_nan() || !(0.0..=1.0).contains(&p) {
             Err(StatsError::BadParams)
         } else {
             Ok(Binomial { p, n })
@@ -279,7 +279,7 @@ impl Discrete<u64, f64> for Binomial {
                 0.0
             }
         } else {
-            (factorial::ln_binomial(self.n as u64, x as u64)
+            (factorial::ln_binomial(self.n, x)
                 + x as f64 * self.p.ln()
                 + (self.n - x) as f64 * (1.0 - self.p).ln())
             .exp()
@@ -310,7 +310,7 @@ impl Discrete<u64, f64> for Binomial {
                 f64::NEG_INFINITY
             }
         } else {
-            factorial::ln_binomial(self.n as u64, x as u64)
+            factorial::ln_binomial(self.n, x)
                 + x as f64 * self.p.ln()
                 + (self.n - x) as f64 * (1.0 - self.p).ln()
         }
diff --git a/src/distribution/cauchy.rs b/src/distribution/cauchy.rs
index dcd81af5..d66a9747 100644
--- a/src/distribution/cauchy.rs
+++ b/src/distribution/cauchy.rs
@@ -138,7 +138,7 @@ impl Max<f64> for Cauchy {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs
index 0a65a8e6..f72376bb 100644
--- a/src/distribution/chi.rs
+++ b/src/distribution/chi.rs
@@ -143,7 +143,7 @@ impl Max<f64> for Chi {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs
index 5551b55f..ab8dc398 100644
--- a/src/distribution/chi_squared.rs
+++ b/src/distribution/chi_squared.rs
@@ -157,7 +157,7 @@ impl Max<f64> for ChiSquared {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/distribution/empirical.rs b/src/distribution/empirical.rs
index 5804f7c1..21e70522 100644
--- a/src/distribution/empirical.rs
+++ b/src/distribution/empirical.rs
@@ -158,14 +158,14 @@ impl ::rand::distributions::Distribution<f64> for Empirical {
 /// Panics if number of samples is zero
 impl Max<f64> for Empirical {
     fn max(&self) -> f64 {
-        self.data.iter().rev().map(|(key, _)| key.0).next().unwrap()
+        self.data.keys().rev().map(|key| key.0) .next().unwrap()
     }
 }
 
 /// Panics if number of samples is zero
 impl Min<f64> for Empirical {
     fn min(&self) -> f64 {
-        self.data.iter().map(|(key, _)| key.0).next().unwrap()
+        self.data.keys().map(|key| key.0).next().unwrap()
     }
 }
 
diff --git a/src/distribution/erlang.rs b/src/distribution/erlang.rs
index e07dff6b..e0721f24 100644
--- a/src/distribution/erlang.rs
+++ b/src/distribution/erlang.rs
@@ -141,7 +141,7 @@ impl Max<f64> for Erlang {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         self.g.max()
diff --git a/src/distribution/exponential.rs b/src/distribution/exponential.rs
index f374b9b2..0856989c 100644
--- a/src/distribution/exponential.rs
+++ b/src/distribution/exponential.rs
@@ -132,7 +132,7 @@ impl Max<f64> for Exp {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/distribution/fisher_snedecor.rs b/src/distribution/fisher_snedecor.rs
index ee70c30e..da8d9570 100644
--- a/src/distribution/fisher_snedecor.rs
+++ b/src/distribution/fisher_snedecor.rs
@@ -174,7 +174,7 @@ impl Max<f64> for FisherSnedecor {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs
index fd993ba5..b8e7df44 100644
--- a/src/distribution/gamma.rs
+++ b/src/distribution/gamma.rs
@@ -2,7 +2,6 @@ use crate::distribution::{Continuous, ContinuousCDF};
 use crate::function::gamma;
 use crate::statistics::*;
 use crate::{Result, StatsError};
-use core::f64::INFINITY as INF;
 use rand::Rng;
 
 /// Implements the [Gamma](https://en.wikipedia.org/wiki/Gamma_distribution)
@@ -172,10 +171,10 @@ impl Max<f64> for Gamma {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
-        INF
+        f64::INFINITY
     }
 }
 
@@ -255,8 +254,8 @@ impl Continuous<f64, f64> for Gamma {
     ///
     /// # Remarks
     ///
-    /// Returns `NAN` if any of `shape` or `rate` are `INF`
-    /// or if `x` is `INF`
+    /// Returns `NAN` if any of `shape` or `rate` are `f64::INFINITY`
+    /// or if `x` is `f64::INFINITY`
     ///
     /// # Formula
     ///
@@ -286,8 +285,8 @@ impl Continuous<f64, f64> for Gamma {
     ///
     /// # Remarks
     ///
-    /// Returns `NAN` if any of `shape` or `rate` are `INF`
-    /// or if `x` is `INF`
+    /// Returns `NAN` if any of `shape` or `rate` are `f64::INFINITY`
+    /// or if `x` is `f64::INFINITY`
     ///
     /// # Formula
     ///
@@ -342,8 +341,8 @@ pub fn sample_unchecked<R: Rng + ?Sized>(rng: &mut R, shape: f64, rate: f64) ->
             };
         }
 
-        v *= v * v;
-        x *= x;
+        v = v * v * v;
+        x = x * x;
         let u: f64 = rng.gen();
         if u < 1.0 - 0.0331 * x * x || u.ln() < 0.5 * x + d * (1.0 - v + v.ln()) {
             return afix * d * v / rate;
@@ -367,7 +366,7 @@ mod tests {
             (1.0, 1.0),
             (10.0, 10.0),
             (10.0, 1.0),
-            (10.0, INF),
+            (10.0, f64::INFINITY),
         ];
 
         for &arg in valid.iter() {
@@ -398,7 +397,7 @@ mod tests {
             ((1.0, 1.0), 1.0),
             ((10.0, 10.0), 1.0),
             ((10.0, 1.0), 10.0),
-            ((10.0, INF), 0.0),
+            ((10.0, f64::INFINITY), 0.0),
         ];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
@@ -413,7 +412,7 @@ mod tests {
             ((1.0, 1.0), 1.0),
             ((10.0, 10.0), 0.1),
             ((10.0, 1.0), 10.0),
-            ((10.0, INF), 0.0),
+            ((10.0, f64::INFINITY), 0.0),
         ];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
@@ -428,7 +427,7 @@ mod tests {
             ((1.0, 1.0), 1.0),
             ((10.0, 10.0), 0.2334690854869339583626209),
             ((10.0, 1.0), 2.53605417848097964238061239),
-            ((10.0, INF), f64::NEG_INFINITY),
+            ((10.0, f64::INFINITY), f64::NEG_INFINITY),
         ];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
@@ -443,7 +442,7 @@ mod tests {
             ((1.0, 1.0), 2.0),
             ((10.0, 10.0), 0.6324555320336758663997787),
             ((10.0, 1.0), 0.63245553203367586639977870),
-            ((10.0, INF), 0.6324555320336758),
+            ((10.0, f64::INFINITY), 0.6324555320336758),
         ];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
@@ -457,7 +456,7 @@ mod tests {
         for &(arg, res) in test.iter() {
             test_case_special(arg, res, 10e-6, f);
         }
-        let test = [((10.0, 10.0), 0.9), ((10.0, 1.0), 9.0), ((10.0, INF), 0.0)];
+        let test = [((10.0, 10.0), 0.9), ((10.0, 1.0), 9.0), ((10.0, f64::INFINITY), 0.0)];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
         }
@@ -471,18 +470,18 @@ mod tests {
             ((1.0, 1.0), 0.0),
             ((10.0, 10.0), 0.0),
             ((10.0, 1.0), 0.0),
-            ((10.0, INF), 0.0),
+            ((10.0, f64::INFINITY), 0.0),
         ];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
         }
         let f = |x: Gamma| x.max();
         let test = [
-            ((1.0, 0.1), INF),
-            ((1.0, 1.0), INF),
-            ((10.0, 10.0), INF),
-            ((10.0, 1.0), INF),
-            ((10.0, INF), INF),
+            ((1.0, 0.1), f64::INFINITY),
+            ((1.0, 1.0), f64::INFINITY),
+            ((10.0, 10.0), f64::INFINITY),
+            ((10.0, 1.0), f64::INFINITY),
+            ((10.0, f64::INFINITY), f64::INFINITY),
         ];
         for &(arg, res) in test.iter() {
             test_case(arg, res, f);
@@ -506,9 +505,9 @@ mod tests {
             test_case(arg, res, f(x));
         }
         //TODO: test special
-        // test_is_nan((10.0, INF), pdf(1.0)); // is this really the behavior we want?
+        // test_is_nan((10.0, f64::INFINITY), pdf(1.0)); // is this really the behavior we want?
         //TODO: test special
-        // (10.0, INF, INF, 0.0, pdf(INF)),];
+        // (10.0, f64::INFINITY, f64::INFINITY, 0.0, pdf(f64::INFINITY)),];
     }
 
     #[test]
@@ -529,13 +528,13 @@ mod tests {
             ((10.0, 10.0), 10.0, -69.0527107131946016148658),
             ((10.0, 1.0), 1.0, -13.8018274800814696112077),
             ((10.0, 1.0), 10.0, -2.07856164313505845504579),
-            ((10.0, INF), INF, f64::NEG_INFINITY),
+            ((10.0, f64::INFINITY), f64::INFINITY, f64::NEG_INFINITY),
         ];
         for &(arg, x, res) in test.iter() {
             test_case(arg, res, f(x));
         }
         // TODO: test special
-        // test_is_nan((10.0, INF), f(1.0)); // is this really the behavior we want?
+        // test_is_nan((10.0, f64::INFINITY), f(1.0)); // is this really the behavior we want?
     }
 
     #[test]
@@ -550,8 +549,8 @@ mod tests {
             ((10.0, 10.0), 10.0, 0.999999999999999999999999),
             ((10.0, 1.0), 1.0, 0.000000111425478338720677),
             ((10.0, 1.0), 10.0, 0.542070285528147791685835),
-            ((10.0, INF), 1.0, 0.0),
-            ((10.0, INF), 10.0, 1.0),
+            ((10.0, f64::INFINITY), 1.0, 0.0),
+            ((10.0, f64::INFINITY), 10.0, 1.0),
         ];
         for &(arg, x, res) in test.iter() {
             test_case(arg, res, f(x));
@@ -575,8 +574,8 @@ mod tests {
             ((10.0, 10.0), 10.0, 1.1253473960842808e-31),
             ((10.0, 1.0), 1.0, 0.9999998885745217),
             ((10.0, 1.0), 10.0, 0.4579297144718528),
-            ((10.0, INF), 1.0, 1.0),
-            ((10.0, INF), 10.0, 0.0),
+            ((10.0, f64::INFINITY), 1.0, 1.0),
+            ((10.0, f64::INFINITY), 10.0, 0.0),
         ];
         for &(arg, x, res) in test.iter() {
             test_case(arg, res, f(x));
diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs
index c7e801c1..f87f5ee0 100644
--- a/src/distribution/geometric.rs
+++ b/src/distribution/geometric.rs
@@ -3,7 +3,7 @@ use crate::statistics::*;
 use crate::{Result, StatsError};
 use rand::distributions::OpenClosed01;
 use rand::Rng;
-use std::{f64, u64};
+use std::f64;
 
 /// Implements the
 /// [Geometric](https://en.wikipedia.org/wiki/Geometric_distribution)
diff --git a/src/distribution/inverse_gamma.rs b/src/distribution/inverse_gamma.rs
index e439be45..31b1d4f6 100644
--- a/src/distribution/inverse_gamma.rs
+++ b/src/distribution/inverse_gamma.rs
@@ -162,7 +162,7 @@ impl Max<f64> for InverseGamma {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/distribution/laplace.rs b/src/distribution/laplace.rs
index 4f73b278..66893b46 100644
--- a/src/distribution/laplace.rs
+++ b/src/distribution/laplace.rs
@@ -173,7 +173,7 @@ impl Max<f64> for Laplace {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
@@ -291,7 +291,6 @@ impl Continuous<f64, f64> for Laplace {
 #[cfg(all(test, feature = "nightly"))]
 mod tests {
     use super::*;
-    use core::f64::INFINITY as INF;
     use rand::thread_rng;
 
     fn try_create(location: f64, scale: f64) -> Laplace {
@@ -349,12 +348,12 @@ mod tests {
     #[test]
     fn test_create() {
         try_create(1.0, 2.0);
-        try_create(-INF, 0.1);
+        try_create(f64::NEG_INFINITY, 0.1);
         try_create(-5.0 - 1.0, 1.0);
         try_create(0.0, 5.0);
         try_create(1.0, 7.0);
         try_create(5.0, 10.0);
-        try_create(INF, INF);
+        try_create(f64::INFINITY, f64::INFINITY);
     }
 
     #[test]
@@ -367,71 +366,71 @@ mod tests {
     #[test]
     fn test_mean() {
         let mean = |x: Laplace| x.mean().unwrap();
-        test_case(-INF, 0.1, -INF, mean);
+        test_case(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, mean);
         test_case(-5.0 - 1.0, 1.0, -6.0, mean);
         test_case(0.0, 5.0, 0.0, mean);
         test_case(1.0, 10.0, 1.0, mean);
-        test_case(INF, INF, INF, mean);
+        test_case(f64::INFINITY, f64::INFINITY, f64::INFINITY, mean);
     }
 
     #[test]
     fn test_variance() {
         let variance = |x: Laplace| x.variance().unwrap();
-        test_almost(-INF, 0.1, 0.02, 1E-12, variance);
+        test_almost(f64::NEG_INFINITY, 0.1, 0.02, 1E-12, variance);
         test_almost(-5.0 - 1.0, 1.0, 2.0, 1E-12, variance);
         test_almost(0.0, 5.0, 50.0, 1E-12, variance);
         test_almost(1.0, 7.0, 98.0, 1E-12, variance);
         test_almost(5.0, 10.0, 200.0, 1E-12, variance);
-        test_almost(INF, INF, INF, 1E-12, variance);
+        test_almost(f64::INFINITY, f64::INFINITY, f64::INFINITY, 1E-12, variance);
     }
     #[test]
     fn test_entropy() {
         let entropy = |x: Laplace| x.entropy().unwrap();
-        test_almost(-INF, 0.1, (2.0 * f64::consts::E * 0.1).ln(), 1E-12, entropy);
+        test_almost(f64::NEG_INFINITY, 0.1, (2.0 * f64::consts::E * 0.1).ln(), 1E-12, entropy);
         test_almost(-6.0, 1.0, (2.0 * f64::consts::E).ln(), 1E-12, entropy);
         test_almost(1.0, 7.0, (2.0 * f64::consts::E * 7.0).ln(), 1E-12, entropy);
         test_almost(5., 10., (2. * f64::consts::E * 10.).ln(), 1E-12, entropy);
-        test_almost(INF, INF, INF, 1E-12, entropy);
+        test_almost(f64::INFINITY, f64::INFINITY, f64::INFINITY, 1E-12, entropy);
     }
 
     #[test]
     fn test_skewness() {
         let skewness = |x: Laplace| x.skewness().unwrap();
-        test_case(-INF, 0.1, 0.0, skewness);
+        test_case(f64::NEG_INFINITY, 0.1, 0.0, skewness);
         test_case(-6.0, 1.0, 0.0, skewness);
         test_case(1.0, 7.0, 0.0, skewness);
         test_case(5.0, 10.0, 0.0, skewness);
-        test_case(INF, INF, 0.0, skewness);
+        test_case(f64::INFINITY, f64::INFINITY, 0.0, skewness);
     }
 
     #[test]
     fn test_mode() {
         let mode = |x: Laplace| x.mode().unwrap();
-        test_case(-INF, 0.1, -INF, mode);
+        test_case(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, mode);
         test_case(-6.0, 1.0, -6.0, mode);
         test_case(1.0, 7.0, 1.0, mode);
         test_case(5.0, 10.0, 5.0, mode);
-        test_case(INF, INF, INF, mode);
+        test_case(f64::INFINITY, f64::INFINITY, f64::INFINITY, mode);
     }
 
     #[test]
     fn test_median() {
         let median = |x: Laplace| x.median();
-        test_case(-INF, 0.1, -INF, median);
+        test_case(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, median);
         test_case(-6.0, 1.0, -6.0, median);
         test_case(1.0, 7.0, 1.0, median);
         test_case(5.0, 10.0, 5.0, median);
-        test_case(INF, INF, INF, median);
+        test_case(f64::INFINITY, f64::INFINITY, f64::INFINITY, median);
     }
 
     #[test]
     fn test_min() {
-        test_case(0.0, 1.0, -INF, |l| l.min());
+        test_case(0.0, 1.0, f64::NEG_INFINITY, |l| l.min());
     }
 
     #[test]
     fn test_max() {
-        test_case(0.0, 1.0, INF, |l| l.max());
+        test_case(0.0, 1.0, f64::INFINITY, |l| l.max());
     }
 
     #[test]
@@ -442,22 +441,22 @@ mod tests {
         test_almost(-1.0, 0.1, 3.8905661205668983e-19, 1E-12, pdf(-5.4));
         test_almost(5.0, 0.1, 5.056107463052243e-43, 1E-12, pdf(-4.9));
         test_almost(-5.0, 0.1, 1.9877248679543235e-30, 1E-12, pdf(2.0));
-        test_almost(INF, 0.1, 0.0, 1E-12, pdf(5.5));
-        test_almost(-INF, 0.1, 0.0, 1E-12, pdf(-0.0));
-        test_almost(0.0, 1.0, 0.0, 1E-12, pdf(INF));
+        test_almost(f64::INFINITY, 0.1, 0.0, 1E-12, pdf(5.5));
+        test_almost(f64::NEG_INFINITY, 0.1, 0.0, 1E-12, pdf(-0.0));
+        test_almost(0.0, 1.0, 0.0, 1E-12, pdf(f64::INFINITY));
         test_almost(1.0, 1.0, 0.00915781944436709, 1E-12, pdf(5.0));
         test_almost(-1.0, 1.0, 0.5, 1E-12, pdf(-1.0));
         test_almost(5.0, 1.0, 0.0012393760883331792, 1E-12, pdf(-1.0));
         test_almost(-5.0, 1.0, 0.0002765421850739168, 1E-12, pdf(2.5));
-        test_almost(INF, 0.1, 0.0, 1E-12, pdf(2.0));
-        test_almost(-INF, 0.1, 0.0, 1E-12, pdf(15.0));
-        test_almost(0.0, INF, 0.0, 1E-12, pdf(89.3));
-        test_almost(1.0, INF, 0.0, 1E-12, pdf(-0.1));
-        test_almost(-1.0, INF, 0.0, 1E-12, pdf(0.1));
-        test_almost(5.0, INF, 0.0, 1E-12, pdf(-6.1));
-        test_almost(-5.0, INF, 0.0, 1E-12, pdf(-10.0));
-        test_is_nan(INF, INF, pdf(2.0));
-        test_is_nan(-INF, INF, pdf(-5.1));
+        test_almost(f64::INFINITY, 0.1, 0.0, 1E-12, pdf(2.0));
+        test_almost(f64::NEG_INFINITY, 0.1, 0.0, 1E-12, pdf(15.0));
+        test_almost(0.0, f64::INFINITY, 0.0, 1E-12, pdf(89.3));
+        test_almost(1.0, f64::INFINITY, 0.0, 1E-12, pdf(-0.1));
+        test_almost(-1.0, f64::INFINITY, 0.0, 1E-12, pdf(0.1));
+        test_almost(5.0, f64::INFINITY, 0.0, 1E-12, pdf(-6.1));
+        test_almost(-5.0, f64::INFINITY, 0.0, 1E-12, pdf(-10.0));
+        test_is_nan(f64::INFINITY, f64::INFINITY, pdf(2.0));
+        test_is_nan(f64::NEG_INFINITY, f64::INFINITY, pdf(-5.1));
     }
 
     #[test]
@@ -468,22 +467,22 @@ mod tests {
         test_almost(-1.0, 0.1, -42.39056208756591, 1E-12, ln_pdf(-5.4));
         test_almost(5.0, 0.1, -97.3905620875659, 1E-12, ln_pdf(-4.9));
         test_almost(-5.0, 0.1, -68.3905620875659, 1E-12, ln_pdf(2.0));
-        test_case(INF, 0.1, -INF, ln_pdf(5.5));
-        test_case(-INF, 0.1, -INF, ln_pdf(-0.0));
-        test_case(0.0, 1.0, -INF, ln_pdf(INF));
+        test_case(f64::INFINITY, 0.1, f64::NEG_INFINITY, ln_pdf(5.5));
+        test_case(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, ln_pdf(-0.0));
+        test_case(0.0, 1.0, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
         test_almost(1.0, 1.0, -4.693147180559945, 1E-12, ln_pdf(5.0));
         test_almost(-1.0, 1.0, -f64::consts::LN_2, 1E-12, ln_pdf(-1.0));
         test_almost(5.0, 1.0, -6.693147180559945, 1E-12, ln_pdf(-1.0));
         test_almost(-5.0, 1.0, -8.193147180559945, 1E-12, ln_pdf(2.5));
-        test_case(INF, 0.1, -INF, ln_pdf(2.0));
-        test_case(-INF, 0.1, -INF, ln_pdf(15.0));
-        test_case(0.0, INF, -INF, ln_pdf(89.3));
-        test_case(1.0, INF, -INF, ln_pdf(-0.1));
-        test_case(-1.0, INF, -INF, ln_pdf(0.1));
-        test_case(5.0, INF, -INF, ln_pdf(-6.1));
-        test_case(-5.0, INF, -INF, ln_pdf(-10.0));
-        test_is_nan(INF, INF, ln_pdf(2.0));
-        test_is_nan(-INF, INF, ln_pdf(-5.1));
+        test_case(f64::INFINITY, 0.1, f64::NEG_INFINITY, ln_pdf(2.0));
+        test_case(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, ln_pdf(15.0));
+        test_case(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(89.3));
+        test_case(1.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-0.1));
+        test_case(-1.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(0.1));
+        test_case(5.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-6.1));
+        test_case(-5.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-10.0));
+        test_is_nan(f64::INFINITY, f64::INFINITY, ln_pdf(2.0));
+        test_is_nan(f64::NEG_INFINITY, f64::INFINITY, ln_pdf(-5.1));
     }
 
     #[test]
diff --git a/src/distribution/log_normal.rs b/src/distribution/log_normal.rs
index 6698dd2a..028869e0 100644
--- a/src/distribution/log_normal.rs
+++ b/src/distribution/log_normal.rs
@@ -136,7 +136,7 @@ impl Max<f64> for LogNormal {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/distribution/negative_binomial.rs b/src/distribution/negative_binomial.rs
index 4c23079c..4c69a869 100644
--- a/src/distribution/negative_binomial.rs
+++ b/src/distribution/negative_binomial.rs
@@ -65,7 +65,7 @@ impl NegativeBinomial {
     /// assert!(result.is_err());
     /// ```
     pub fn new(r: f64, p: f64) -> Result<NegativeBinomial> {
-        if p.is_nan() || p < 0.0 || p > 1.0 || r.is_nan() || r < 0.0 {
+        if p.is_nan() || !(0.0..=1.0).contains(&p) || r.is_nan() || r < 0.0 {
             Err(StatsError::BadParams)
         } else {
             Ok(NegativeBinomial { r, p })
@@ -173,7 +173,7 @@ impl Max<u64> for NegativeBinomial {
     /// u64::MAX
     /// ```
     fn max(&self) -> u64 {
-        std::u64::MAX
+        u64::MAX
     }
 }
 
diff --git a/src/distribution/normal.rs b/src/distribution/normal.rs
index d1cc91f2..bcd01b81 100644
--- a/src/distribution/normal.rs
+++ b/src/distribution/normal.rs
@@ -130,7 +130,7 @@ impl Min<f64> for Normal {
     /// # Formula
     ///
     /// ```text
-    /// -INF
+    /// f64::NEG_INFINITY
     /// ```
     fn min(&self) -> f64 {
         f64::NEG_INFINITY
@@ -144,7 +144,7 @@ impl Max<f64> for Normal {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/distribution/pareto.rs b/src/distribution/pareto.rs
index eca17f82..031205eb 100644
--- a/src/distribution/pareto.rs
+++ b/src/distribution/pareto.rs
@@ -160,7 +160,7 @@ impl Max<f64> for Pareto {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
@@ -174,7 +174,7 @@ impl Distribution<f64> for Pareto {
     ///
     /// ```text
     /// if α <= 1 {
-    ///     INF
+    ///     f64::INFINITY
     /// } else {
     ///     (α * x_m)/(α - 1)
     /// }
@@ -194,7 +194,7 @@ impl Distribution<f64> for Pareto {
     ///
     /// ```text
     /// if α <= 2 {
-    ///     INF
+    ///     f64::INFINITY
     /// } else {
     ///     (x_m/(α - 1))^2 * (α/(α - 2))
     /// }
@@ -308,7 +308,7 @@ impl Continuous<f64, f64> for Pareto {
     ///
     /// ```text
     /// if x < x_m {
-    ///     -INF
+    ///     f64::NEG_INFINITY
     /// } else {
     ///     ln(α) + α*ln(x_m) - (α + 1)*ln(x)
     /// }
diff --git a/src/distribution/poisson.rs b/src/distribution/poisson.rs
index 47d614ce..e8f98bef 100644
--- a/src/distribution/poisson.rs
+++ b/src/distribution/poisson.rs
@@ -4,7 +4,6 @@ use crate::statistics::*;
 use crate::{Result, StatsError};
 use rand::Rng;
 use std::f64;
-use std::u64;
 
 /// Implements the [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution)
 /// distribution
@@ -234,7 +233,7 @@ impl Discrete<u64, f64> for Poisson {
     ///
     /// where `λ` is the rate
     fn pmf(&self, x: u64) -> f64 {
-        (-self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x as u64)).exp()
+        (-self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x)).exp()
     }
 
     /// Calculates the log probability mass function for the poisson
@@ -249,7 +248,7 @@ impl Discrete<u64, f64> for Poisson {
     ///
     /// where `λ` is the rate
     fn ln_pmf(&self, x: u64) -> f64 {
-        -self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x as u64)
+        -self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x)
     }
 }
 /// Generates one sample from the Poisson distribution either by
diff --git a/src/distribution/students_t.rs b/src/distribution/students_t.rs
index 21981513..4f84c489 100644
--- a/src/distribution/students_t.rs
+++ b/src/distribution/students_t.rs
@@ -210,7 +210,7 @@ impl Min<f64> for StudentsT {
     /// # Formula
     ///
     /// ```text
-    /// -INF
+    /// f64::NEG_INFINITY
     /// ```
     fn min(&self) -> f64 {
         f64::NEG_INFINITY
@@ -224,7 +224,7 @@ impl Max<f64> for StudentsT {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
@@ -261,7 +261,7 @@ impl Distribution<f64> for StudentsT {
     /// # Formula
     ///
     /// ```text
-    /// if v == INF {
+    /// if v == f64::INFINITY {
     ///     Some(σ^2)
     /// } else if freedom > 2.0 {
     ///     Some(v * σ^2 / (v - 2))
diff --git a/src/distribution/weibull.rs b/src/distribution/weibull.rs
index c414c20c..4f04403d 100644
--- a/src/distribution/weibull.rs
+++ b/src/distribution/weibull.rs
@@ -156,7 +156,7 @@ impl Max<f64> for Weibull {
     /// # Formula
     ///
     /// ```text
-    /// INF
+    /// f64::INFINITY
     /// ```
     fn max(&self) -> f64 {
         f64::INFINITY
diff --git a/src/function/beta.rs b/src/function/beta.rs
index 23fb3430..fec184f2 100644
--- a/src/function/beta.rs
+++ b/src/function/beta.rs
@@ -327,11 +327,7 @@ pub fn inv_beta_reg(mut a: f64, mut b: f64, mut x: f64) -> f64 {
         }
     }
 
-    if p < 0.0001 {
-        p = 0.0001;
-    } else if 0.9999 < p {
-        p = 0.9999;
-    }
+    p = p.clamp(0.0001, 0.9999);
 
     // Remark AS R83
     // http://www.jstor.org/stable/2347779
@@ -365,7 +361,7 @@ pub fn inv_beta_reg(mut a: f64, mut b: f64, mut x: f64) -> f64 {
 
                 if sq < prev {
                     pnext = p - adj;
-                    if 0.0 <= pnext && pnext <= 1.0 {
+                    if (0.0..=1.0).contains(&pnext) {
                         break;
                     }
                 }
diff --git a/src/function/factorial.rs b/src/function/factorial.rs
index eac59bcf..e06e9bcf 100644
--- a/src/function/factorial.rs
+++ b/src/function/factorial.rs
@@ -4,7 +4,6 @@
 use crate::error::StatsError;
 use crate::function::gamma;
 use crate::Result;
-use core::f64::INFINITY as INF;
 
 /// The maximum factorial representable
 /// by a 64-bit floating point without
@@ -20,7 +19,7 @@ pub const MAX_FACTORIAL: usize = 170;
 /// Returns `f64::INFINITY` if `x > 170`
 pub fn factorial(x: u64) -> f64 {
     let x = x as usize;
-    FCACHE.get(x).map_or(INF, |&fac| fac)
+    FCACHE.get(x).map_or(f64::INFINITY, |&fac| fac)
 }
 
 /// Computes the logarithmic factorial function `x -> ln(x!)`
@@ -124,8 +123,8 @@ mod tests {
 
     #[test]
     fn test_factorial_overflow() {
-        assert_eq!(factorial(172), INF);
-        assert_eq!(factorial(u64::MAX), INF);
+        assert_eq!(factorial(172), f64::INFINITY);
+        assert_eq!(factorial(u64::MAX), f64::INFINITY);
     }
 
     #[test]
diff --git a/src/generate.rs b/src/generate.rs
index 1f6102d2..59783598 100644
--- a/src/generate.rs
+++ b/src/generate.rs
@@ -323,7 +323,7 @@ impl InfiniteSawtooth {
                 0.0,
                 delay,
             ),
-            low_value: low_value as f64,
+            low_value,
         }
     }
 }