diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs
index 208047db..ead08266 100644
--- a/src/distribution/chi.rs
+++ b/src/distribution/chi.rs
@@ -119,12 +119,13 @@ impl ContinuousCDF<f64, f64> for Chi {
     /// where `k` is the degrees of freedom and `P` is
     /// the regularized lower incomplete Gamma function
     fn cdf(&self, x: f64) -> f64 {
-        if x == f64::INFINITY {
-            1.0
-        } else if x <= 0.0 {
+        let u = x * x / 2.0;
+        if x <= 0.0 || u == 0.0 {
             0.0
+        } else if u == f64::INFINITY {
+            1.0
         } else {
-            gamma::gamma_lr(self.freedom() as f64 / 2.0, x * x / 2.0)
+            gamma::gamma_lr(self.freedom() as f64 / 2.0, u)
         }
     }
 
@@ -140,12 +141,13 @@ impl ContinuousCDF<f64, f64> for Chi {
     /// where `k` is the degrees of freedom and `P` is
     /// the regularized upper incomplete Gamma function
     fn sf(&self, x: f64) -> f64 {
-        if x == f64::INFINITY {
-            0.0
-        } else if x <= 0.0 {
+        let u = x * x / 2.0;
+        if x <= 0.0 || u == 0.0 {
             1.0
+        } else if u == f64::INFINITY {
+            0.0
         } else {
-            gamma::gamma_ur(self.freedom() as f64 / 2.0, x * x / 2.0)
+            gamma::gamma_ur(self.freedom() as f64 / 2.0, u)
         }
     }
 
@@ -503,6 +505,11 @@ mod tests {
         test_exact(2, 1.0, cdf(f64::INFINITY));
         test_exact(2, 0.0, cdf(0.0));
         test_exact(2, 1.0, cdf(f64::INFINITY));
+        // x * x / 2 underflows to 0 or overflows to infinity
+        test_exact(1, 0.0, cdf(1e-200));
+        test_exact(3, 0.0, cdf(1e-200));
+        test_exact(1, 1.0, cdf(1e160));
+        test_exact(3, 1.0, cdf(f64::MAX));
     }
 
     #[test]
@@ -520,18 +527,23 @@ mod tests {
         test_exact(2, 0.0, sf(f64::INFINITY));
         test_exact(2, 1.0, sf(0.0));
         test_exact(2, 0.0, sf(f64::INFINITY));
+        test_exact(1, 1.0, sf(1e-200));
+        test_exact(1, 0.0, sf(1e160));
+        test_exact(3, 0.0, sf(f64::MAX));
     }
 
     #[test]
     fn test_neg_cdf() {
         let cdf = |arg: f64| move |x: Chi| x.cdf(arg);
         test_exact(1, 0.0, cdf(-1.0));
+        test_exact(1, 0.0, cdf(f64::NEG_INFINITY));
     }
 
     #[test]
     fn test_neg_sf() {
         let sf = |arg: f64| move |x: Chi| x.sf(arg);
         test_exact(1, 1.0, sf(-1.0));
+        test_exact(1, 1.0, sf(f64::NEG_INFINITY));
     }
 
     #[test]
diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs
index 6df26e38..a4d2da93 100644
--- a/src/distribution/gamma.rs
+++ b/src/distribution/gamma.rs
@@ -152,8 +152,10 @@ impl ContinuousCDF<f64, f64> for Gamma {
             1.0
         } else if self.rate.is_infinite() {
             0.0
-        } else if x.is_infinite() {
+        } else if x.is_infinite() || (x * self.rate).is_infinite() {
             1.0
+        } else if x * self.rate == 0.0 {
+            0.0
         } else {
             gamma::gamma_lr(self.shape, x * self.rate)
         }
@@ -177,8 +179,10 @@ impl ContinuousCDF<f64, f64> for Gamma {
             0.0
         } else if self.rate.is_infinite() {
             1.0
-        } else if x.is_infinite() {
+        } else if x.is_infinite() || (x * self.rate).is_infinite() {
             0.0
+        } else if x * self.rate == 0.0 {
+            1.0
         } else {
             gamma::gamma_ur(self.shape, x * self.rate)
         }
@@ -745,6 +749,14 @@ mod tests {
         test_relative(1.0, 0.1, 0.0, |x| x.cdf(0.0));
     }
 
+    #[test]
+    fn test_cdf_sf_when_x_times_rate_overflows_or_underflows() {
+        test_exact(2.5, 1.5, 1.0, |g| g.cdf(f64::MAX));
+        test_exact(2.5, 1.5, 0.0, |g| g.sf(f64::MAX));
+        test_exact(1.0, 0.5, 0.0, |g| g.cdf(5e-324));
+        test_exact(1.0, 0.5, 1.0, |g| g.sf(5e-324));
+    }
+
     #[test]
     fn test_infinite_rate_special_point_is_exact() {
         let dist = create_ok(10.0, f64::INFINITY);