diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs
index 38ae1517..c8e9a0f5 100644
--- a/src/distribution/beta.rs
+++ b/src/distribution/beta.rs
@@ -1,4 +1,4 @@
-use crate::distribution::{Continuous, ContinuousCDF, InverseCdfError};
+use crate::distribution::{Continuous, ContinuousCDF};
 use crate::function::{beta, gamma};
 use crate::statistics::*;
 #[cfg(not(feature = "std"))]
@@ -181,10 +181,6 @@ impl ContinuousCDF<f64, f64> for Beta {
     /// Calculates the inverse cumulative distribution function for the beta
     /// distribution at `x`.
     ///
-    /// # Panics
-    ///
-    /// If x is not in `[0, 1]`.
-    ///
     /// # Formula
     ///
     /// ```text
@@ -194,34 +190,7 @@ impl ContinuousCDF<f64, f64> for Beta {
     /// where `α` is shapeA, `β` is shapeB, and `I_x` is the inverse of the
     /// regularized lower incomplete beta function.
     fn inverse_cdf(&self, x: f64) -> f64 {
-        if !(0.0..=1.0).contains(&x) {
-            panic!("x must be in [0, 1]");
-        } else {
-            beta::inv_beta_reg(self.shape_a, self.shape_b, x)
-        }
-    }
-
-    /// Calculates the inverse cumulative distribution function for the beta
-    /// distribution at `x`.
-    ///
-    /// # Returns an error instead of a panic
-    ///
-    /// If x is not in `[0, 1]`.
-    ///
-    /// # Formula
-    ///
-    /// ```text
-    /// I^{-1}_x(α, β)
-    /// ```
-    ///
-    /// where `α` is shapeA, `β` is shapeB, and `I_x` is the inverse of the
-    /// regularized lower incomplete beta function.
-    fn try_inverse_cdf(&self, x: f64) -> Result<f64, InverseCdfError> {
-        if !(0.0..=1.0).contains(&x) {
-            Err(InverseCdfError::ArgumentOutOfRange)
-        } else {
-            Ok(beta::inv_beta_reg(self.shape_a, self.shape_b, x))
-        }
+        beta::inv_beta_reg(self.shape_a, self.shape_b, x)
     }
 }
 
diff --git a/src/distribution/categorical.rs b/src/distribution/categorical.rs
index c4de36f0..a0187cf4 100644
--- a/src/distribution/categorical.rs
+++ b/src/distribution/categorical.rs
@@ -1,4 +1,4 @@
-use crate::distribution::{Discrete, DiscreteCDF};
+use crate::distribution::{Discrete, DiscreteCDF, DiscreteInverseCdfError};
 use crate::statistics::*;
 use alloc::vec::Vec;
 #[cfg(not(feature = "std"))]
@@ -181,10 +181,6 @@ impl DiscreteCDF<u64, f64> for Categorical {
     /// categorical
     /// distribution at `x`
     ///
-    /// # Panics
-    ///
-    /// If `x <= 0.0` or `x >= 1.0`
-    ///
     /// # Formula
     ///
     /// ```text
@@ -195,12 +191,23 @@ impl DiscreteCDF<u64, f64> for Categorical {
     /// and `f(x)` is defined as `p_x + f(x - 1)` and `f(0) = p_0` where
     /// `p_x` is the `x`th probability mass
     fn inverse_cdf(&self, x: f64) -> u64 {
-        if x >= 1.0 || x <= 0.0 {
-            panic!("x must be in [0, 1]")
-        }
-
         self.locate(x)
     }
+
+    /// Calculates the inverse cumulative distribution function for the
+    /// categorical distribution at `x`, returning an error instead of
+    /// panicking if `x` is not in the open interval `(0.0, 1.0)`, including
+    /// when `x` is NaN.
+    fn try_inverse_cdf(&self, x: f64) -> Result<u64, DiscreteInverseCdfError<f64>> {
+        // Written as a conjunction of `>`/`<`, not `x >= 1.0 || x <= 0.0`, so
+        // that NaN (which compares false against everything) is rejected
+        // rather than falling through to `locate`, which panics on NaN.
+        if x > 0.0 && x < 1.0 {
+            Ok(self.inverse_cdf(x))
+        } else {
+            Err(DiscreteInverseCdfError::ArgumentOutOfRange)
+        }
+    }
 }
 
 impl Min<u64> for Categorical {
@@ -517,17 +524,30 @@ mod tests {
     }
 
     #[test]
-    #[should_panic]
-    fn test_inverse_cdf_input_low() {
+    fn test_try_inverse_cdf_input_low() {
+        let dist = create_ok(&[4.0, 2.5, 2.5, 1.0]);
+        assert_eq!(
+            dist.try_inverse_cdf(0.0),
+            Err(DiscreteInverseCdfError::ArgumentOutOfRange)
+        );
+    }
+
+    #[test]
+    fn test_try_inverse_cdf_input_high() {
         let dist = create_ok(&[4.0, 2.5, 2.5, 1.0]);
-        dist.inverse_cdf(0.0);
+        assert_eq!(
+            dist.try_inverse_cdf(1.0),
+            Err(DiscreteInverseCdfError::ArgumentOutOfRange)
+        );
     }
 
     #[test]
-    #[should_panic]
-    fn test_inverse_cdf_input_high() {
+    fn test_try_inverse_cdf_nan_does_not_panic() {
         let dist = create_ok(&[4.0, 2.5, 2.5, 1.0]);
-        dist.inverse_cdf(1.0);
+        assert_eq!(
+            dist.try_inverse_cdf(f64::NAN),
+            Err(DiscreteInverseCdfError::ArgumentOutOfRange)
+        );
     }
 
     #[test]
diff --git a/src/distribution/cauchy.rs b/src/distribution/cauchy.rs
index 0d978b6f..f0ce40ec 100644
--- a/src/distribution/cauchy.rs
+++ b/src/distribution/cauchy.rs
@@ -163,11 +163,7 @@ impl ContinuousCDF<f64, f64> for Cauchy {
     ///
     /// where `x_0` is the location and `γ` is the scale
     fn inverse_cdf(&self, x: f64) -> f64 {
-        if !(0.0..=1.0).contains(&x) {
-            panic!("x must be in [0, 1]");
-        } else {
-            self.location + self.scale * (f64_consts::PI * (x - 0.5)).tan()
-        }
+        self.location + self.scale * (f64_consts::PI * (x - 0.5)).tan()
     }
 }
 
diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs
index 208047db..c69ef696 100644
--- a/src/distribution/chi.rs
+++ b/src/distribution/chi.rs
@@ -151,14 +151,7 @@ impl ContinuousCDF<f64, f64> for Chi {
 
     /// Calculates the inverse cumulative distribution function for the chi
     /// distribution at `p`, i.e. the `p`-quantile.
-    ///
-    /// # Panics
-    ///
-    /// If `p` is not in `[0, 1]`.
     fn inverse_cdf(&self, p: f64) -> f64 {
-        if !(0.0..=1.0).contains(&p) {
-            panic!("p must be in [0, 1]")
-        }
         if p == 0.0 {
             return self.min();
         }
@@ -367,6 +360,7 @@ impl Continuous<f64, f64> for Chi {
 mod tests {
     use super::*;
     use crate::distribution::internal::density_util;
+    use crate::distribution::InverseCdfError;
     crate::distribution::internal::testing_boiler!(freedom: u64; Chi; ChiError);
 
     #[test]
@@ -588,14 +582,18 @@ mod tests {
     }
 
     #[test]
-    #[should_panic(expected = "p must be in [0, 1]")]
-    fn test_inverse_cdf_p_above_one() {
-        create_ok(3).inverse_cdf(1.0 + f64::EPSILON);
+    fn test_try_inverse_cdf_p_above_one() {
+        assert_eq!(
+            create_ok(3).try_inverse_cdf(1.0 + f64::EPSILON),
+            Err(InverseCdfError::ArgumentOutOfRange)
+        );
     }
 
     #[test]
-    #[should_panic(expected = "p must be in [0, 1]")]
-    fn test_inverse_cdf_p_below_zero() {
-        create_ok(3).inverse_cdf(-1e-300);
+    fn test_try_inverse_cdf_p_below_zero() {
+        assert_eq!(
+            create_ok(3).try_inverse_cdf(-1e-300),
+            Err(InverseCdfError::ArgumentOutOfRange)
+        );
     }
 }
diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs
index 882a3772..8c24898d 100644
--- a/src/distribution/chi_squared.rs
+++ b/src/distribution/chi_squared.rs
@@ -1,4 +1,4 @@
-use crate::distribution::{Continuous, ContinuousCDF, Gamma, GammaError};
+use crate::distribution::{Continuous, ContinuousCDF, Gamma, GammaError, InverseCdfError};
 use crate::statistics::*;
 
 /// Implements the
@@ -152,6 +152,10 @@ impl ContinuousCDF<f64, f64> for ChiSquared {
     fn inverse_cdf(&self, p: f64) -> f64 {
         self.g.inverse_cdf(p)
     }
+
+    fn try_inverse_cdf(&self, p: f64) -> Result<f64, InverseCdfError> {
+        self.g.try_inverse_cdf(p)
+    }
 }
 
 impl Min<f64> for ChiSquared {
diff --git a/src/distribution/erlang.rs b/src/distribution/erlang.rs
index 8fd5440f..fa2aa05b 100644
--- a/src/distribution/erlang.rs
+++ b/src/distribution/erlang.rs
@@ -1,4 +1,4 @@
-use crate::distribution::{Continuous, ContinuousCDF, Gamma, GammaError};
+use crate::distribution::{Continuous, ContinuousCDF, Gamma, GammaError, InverseCdfError};
 use crate::statistics::*;
 
 /// Implements the [Erlang](https://en.wikipedia.org/wiki/Erlang_distribution)
@@ -137,6 +137,10 @@ impl ContinuousCDF<f64, f64> for Erlang {
     fn inverse_cdf(&self, p: f64) -> f64 {
         self.g.inverse_cdf(p)
     }
+
+    fn try_inverse_cdf(&self, p: f64) -> Result<f64, InverseCdfError> {
+        self.g.try_inverse_cdf(p)
+    }
 }
 
 impl Min<f64> for Erlang {
diff --git a/src/distribution/fisher_snedecor.rs b/src/distribution/fisher_snedecor.rs
index 334e19b2..7d868f98 100644
--- a/src/distribution/fisher_snedecor.rs
+++ b/src/distribution/fisher_snedecor.rs
@@ -201,12 +201,8 @@ impl ContinuousCDF<f64, f64> for FisherSnedecor {
     /// the second degree of freedom, and `I` is the regularized incomplete
     /// beta function
     fn inverse_cdf(&self, x: f64) -> f64 {
-        if !(0.0..=1.0).contains(&x) {
-            panic!("x must be in [0, 1]");
-        } else {
-            let z = beta::inv_beta_reg(self.freedom_1 / 2.0, self.freedom_2 / 2.0, x);
-            self.freedom_2 / (self.freedom_1 * (1.0 / z - 1.0))
-        }
+        let z = beta::inv_beta_reg(self.freedom_1 / 2.0, self.freedom_2 / 2.0, x);
+        self.freedom_2 / (self.freedom_1 * (1.0 / z - 1.0))
     }
 }
 
diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs
index 6c376774..67a1edf8 100644
--- a/src/distribution/gamma.rs
+++ b/src/distribution/gamma.rs
@@ -185,9 +185,6 @@ impl ContinuousCDF<f64, f64> for Gamma {
     }
 
     fn inverse_cdf(&self, p: f64) -> f64 {
-        if !(0.0..=1.0).contains(&p) {
-            panic!("default inverse_cdf implementation should be provided probability on [0,1]")
-        }
         if p == 0.0 {
             return self.min();
         };
diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs
index 2972081e..6f2eb725 100644
--- a/src/distribution/geometric.rs
+++ b/src/distribution/geometric.rs
@@ -1,4 +1,4 @@
-use crate::distribution::{Discrete, DiscreteCDF};
+use crate::distribution::{Discrete, DiscreteCDF, DiscreteInverseCdfError};
 use crate::statistics::*;
 use core::f64::consts as f64_consts;
 #[cfg(not(feature = "std"))]
@@ -184,43 +184,58 @@ impl DiscreteCDF<u64, f64> for Geometric {
     ///
     /// # Panics
     ///
-    /// Panics if `x` is not in `[0, 1]`.
-    /// Panics if the result would exceed `u64::MAX` (p is too small for the given x).
-    /// Panics if intermediate f64 computation overflows (p is pathologically small).
+    /// Panics if `x` is not in `[0, 1]` (including NaN), if the result would
+    /// exceed `u64::MAX` (p is too small for the given x), or if intermediate
+    /// f64 computation overflows (p is pathologically small).
     fn inverse_cdf(&self, x: f64) -> u64 {
+        match self.try_inverse_cdf(x) {
+            Ok(k) => k,
+            Err(DiscreteInverseCdfError::ArgumentOutOfRange) => {
+                panic!("inverse_cdf: x must be in [0, 1], was {x}")
+            }
+            Err(DiscreteInverseCdfError::NotRepresentable(k)) if !k.is_finite() => panic!(
+                "inverse_cdf: intermediate value overflowed f64; p ({}) is too small",
+                self.p
+            ),
+            Err(DiscreteInverseCdfError::NotRepresentable(_)) => panic!(
+                "inverse_cdf: result exceeds u64::MAX; p ({}) is too small for x ({x})",
+                self.p
+            ),
+        }
+    }
+
+    /// Calculates the inverse cumulative distribution function for the
+    /// geometric distribution at `x`, returning an error instead of panicking
+    /// if `x` is not in `[0, 1]` (including NaN), or if the exact quantile
+    /// is not representable as a `u64`.
+    fn try_inverse_cdf(&self, x: f64) -> Result<u64, DiscreteInverseCdfError<f64>> {
+        if !(0.0..=1.0).contains(&x) {
+            return Err(DiscreteInverseCdfError::ArgumentOutOfRange);
+        }
         if x == <f64>::zero() {
-            return self.min();
+            return Ok(self.min());
         }
         if x == <f64>::one() {
             // iCDF(1) = +∞; saturate to supremum of u64 domain
-            return self.max();
-        }
-        if !(<f64>::zero()..=<f64>::one()).contains(&x) {
-            panic!("x must be in [0, 1]");
+            return Ok(self.max());
         }
         if self.p == 1.0 {
             // degenerate distribution: all mass at k=1
-            return self.min();
+            return Ok(self.min());
         }
         // cdf(1) = p exactly, so every probability in (0, p] maps to the mode.
         // Handle this before the closed form so platform-dependent ln1p/expm1
         // noise in cdf cannot push the answer to 2 (observed on Windows MSVC).
         if x <= self.p {
-            return self.min();
+            return Ok(self.min());
         }
         let k = (-x).ln_1p() / (-self.p).ln_1p();
         if !k.is_finite() {
-            panic!(
-                "inverse_cdf: intermediate value overflowed f64; p ({}) is too small",
-                self.p
-            );
+            return Err(DiscreteInverseCdfError::NotRepresentable(k));
         }
         let k = k.ceil();
         if k >= u64::MAX as f64 {
-            panic!(
-                "inverse_cdf: result exceeds u64::MAX; p ({}) is too small for x ({})",
-                self.p, x
-            );
+            return Err(DiscreteInverseCdfError::NotRepresentable(k));
         }
 
         // `ln1p(-x) / ln1p(-p)` is only approximately integral at the step
@@ -233,13 +248,13 @@ impl DiscreteCDF<u64, f64> for Geometric {
         // Ordinary rounding puts the closed form within one step, so this is
         // the path essentially always taken: two or three cdf evaluations.
         if is_answer(candidate) {
-            return candidate;
+            return Ok(candidate);
         }
         if candidate > self.min() && is_answer(candidate - 1) {
-            return candidate - 1;
+            return Ok(candidate - 1);
         }
         if candidate < u64::MAX && is_answer(candidate + 1) {
-            return candidate + 1;
+            return Ok(candidate + 1);
         }
 
         // Once x is within a few ulp of 1, `1 - x` has lost significant bits and
@@ -266,7 +281,7 @@ impl DiscreteCDF<u64, f64> for Geometric {
                 lo = mid + 1;
             }
         }
-        lo
+        Ok(lo)
     }
 }
 
@@ -421,6 +436,7 @@ impl Discrete<u64, f64> for Geometric {
 mod tests {
     use super::*;
     use crate::distribution::internal::density_util;
+    use crate::distribution::DiscreteInverseCdfError;
     use crate::prec;
 
     crate::distribution::internal::testing_boiler!(p: f64; Geometric; GeometricError);
@@ -748,10 +764,37 @@ mod tests {
     }
 
     #[test]
-    #[should_panic]
-    fn test_inverse_cdf_panic() {
-        let invcdf = |arg: f64| move |x: Geometric| x.inverse_cdf(arg);
-        test_exact(1., 1, invcdf(2.));
+    fn test_try_inverse_cdf_out_of_range() {
+        assert_eq!(
+            create_ok(0.5).try_inverse_cdf(2.),
+            Err(DiscreteInverseCdfError::ArgumentOutOfRange)
+        );
+    }
+
+    #[test]
+    fn test_try_inverse_cdf_nan_does_not_panic() {
+        assert_eq!(
+            create_ok(0.5).try_inverse_cdf(f64::NAN),
+            Err(DiscreteInverseCdfError::ArgumentOutOfRange)
+        );
+    }
+
+    #[test]
+    fn test_try_inverse_cdf_intermediate_overflow_not_representable() {
+        let distribution = Geometric::new(f64::from_bits(1)).unwrap();
+        assert!(matches!(
+            distribution.try_inverse_cdf(0.5),
+            Err(DiscreteInverseCdfError::NotRepresentable(k)) if !k.is_finite()
+        ));
+    }
+
+    #[test]
+    fn test_try_inverse_cdf_result_overflow_not_representable() {
+        let distribution = Geometric::new(1e-20).unwrap();
+        assert!(matches!(
+            distribution.try_inverse_cdf(0.5),
+            Err(DiscreteInverseCdfError::NotRepresentable(k)) if k.is_finite()
+        ));
     }
 
     #[test]
diff --git a/src/distribution/inverse_gamma.rs b/src/distribution/inverse_gamma.rs
index 3630b072..06dd78f4 100644
--- a/src/distribution/inverse_gamma.rs
+++ b/src/distribution/inverse_gamma.rs
@@ -174,14 +174,7 @@ impl ContinuousCDF<f64, f64> for InverseGamma {
 
     /// Calculates the inverse cumulative distribution function for the inverse
     /// gamma distribution at `p`, i.e. the `p`-quantile.
-    ///
-    /// # Panics
-    ///
-    /// If `p` is not in `[0, 1]`.
     fn inverse_cdf(&self, p: f64) -> f64 {
-        if !(0.0..=1.0).contains(&p) {
-            panic!("p must be in [0, 1]")
-        }
         if p == 0.0 {
             return self.min();
         }
@@ -370,6 +363,7 @@ impl Continuous<f64, f64> for InverseGamma {
 mod tests {
     use super::*;
     use crate::distribution::internal::density_util;
+    use crate::distribution::InverseCdfError;
     use crate::prec;
     use core::f64::consts as f64_consts;
 
@@ -587,14 +581,18 @@ mod tests {
     }
 
     #[test]
-    #[should_panic(expected = "p must be in [0, 1]")]
-    fn test_inverse_cdf_p_above_one() {
-        create_ok(1.0, 1.0).inverse_cdf(1.0 + f64::EPSILON);
+    fn test_try_inverse_cdf_p_above_one() {
+        assert_eq!(
+            create_ok(1.0, 1.0).try_inverse_cdf(1.0 + f64::EPSILON),
+            Err(InverseCdfError::ArgumentOutOfRange)
+        );
     }
 
     #[test]
-    #[should_panic(expected = "p must be in [0, 1]")]
-    fn test_inverse_cdf_p_below_zero() {
-        create_ok(1.0, 1.0).inverse_cdf(-1e-300);
+    fn test_try_inverse_cdf_p_below_zero() {
+        assert_eq!(
+            create_ok(1.0, 1.0).try_inverse_cdf(-1e-300),
+            Err(InverseCdfError::ArgumentOutOfRange)
+        );
     }
 }
diff --git a/src/distribution/laplace.rs b/src/distribution/laplace.rs
index 3543dcaf..13edcf5d 100644
--- a/src/distribution/laplace.rs
+++ b/src/distribution/laplace.rs
@@ -167,9 +167,6 @@ impl ContinuousCDF<f64, f64> for Laplace {
     ///
     /// where `μ` is the location, `b` is the scale
     fn inverse_cdf(&self, p: f64) -> f64 {
-        if p <= 0. || 1. <= p {
-            panic!("p must be in [0, 1]");
-        };
         if p <= 0.5 {
             self.location + self.scale * (2. * p).ln()
         } else {
diff --git a/src/distribution/levy.rs b/src/distribution/levy.rs
index 8102bd75..067828b2 100644
--- a/src/distribution/levy.rs
+++ b/src/distribution/levy.rs
@@ -166,10 +166,6 @@ impl ContinuousCDF<f64, f64> for Levy {
     /// Calculates the inverse cumulative distribution function for the
     /// normal distribution at `x`.
     ///
-    /// # Panics
-    ///
-    /// If `x < 0.0` or `x > 1.0`
-    ///
     /// # Formula
     ///
     /// ```text
@@ -179,11 +175,7 @@ impl ContinuousCDF<f64, f64> for Levy {
     /// where `μ` is the mean, `σ` is the standard deviation and `erfc_inv` is
     /// the inverse of the complementary error function
     fn inverse_cdf(&self, x: f64) -> f64 {
-        if !(0.0..=1.0).contains(&x) {
-            panic!("x must be in [0, 1]");
-        } else {
-            self.mu + 0.5 * self.c / (erfc_inv(x).powf(2.0))
-        }
+        self.mu + 0.5 * self.c / (erfc_inv(x).powf(2.0))
     }
 }
 
diff --git a/src/distribution/log_normal.rs b/src/distribution/log_normal.rs
index f8e9e83c..7953362d 100644
--- a/src/distribution/log_normal.rs
+++ b/src/distribution/log_normal.rs
@@ -181,10 +181,6 @@ impl ContinuousCDF<f64, f64> for LogNormal {
     /// Calculates the inverse cumulative distribution function for the
     /// log-normal distribution at `p`
     ///
-    /// # Panics
-    ///
-    /// If `p < 0.0` or `p > 1.0`
-    ///
     /// # Formula
     ///
     /// ```text
@@ -198,10 +194,8 @@ impl ContinuousCDF<f64, f64> for LogNormal {
             0.0
         } else if p < 1.0 {
             (self.location - (self.scale * f64_consts::SQRT_2 * erf::erfc_inv(2.0 * p))).exp()
-        } else if p == 1.0 {
-            f64::INFINITY
         } else {
-            panic!("p must be within [0.0, 1.0]");
+            f64::INFINITY
         }
     }
 }
diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs
index c31b53d7..83c155df 100644
--- a/src/distribution/mod.rs
+++ b/src/distribution/mod.rs
@@ -115,6 +115,39 @@ impl core::fmt::Display for InverseCdfError {
 
 impl core::error::Error for InverseCdfError {}
 
+/// Represents the errors that can occur when computing [`DiscreteCDF::try_inverse_cdf`].
+///
+/// Unlike [`InverseCdfError`], a discrete quantile can fail for a second reason
+/// beyond an out-of-range argument: the exact answer may not be representable in
+/// the distribution's variate type `K` (for example a `u64`), because it is NaN,
+/// non-finite, or a finite value outside `K`'s range. `NotRepresentable` carries
+/// that underlying value for diagnostics.
+#[derive(Copy, Clone, PartialEq, Debug)]
+#[non_exhaustive]
+pub enum DiscreteInverseCdfError<T> {
+    /// The argument `p` is outside the closed interval `[0, 1]`, or is NaN.
+    ArgumentOutOfRange,
+    /// The exact quantile is not representable in the distribution's variate
+    /// type. Carries the floating-point value that could not be converted.
+    NotRepresentable(T),
+}
+
+impl<T: core::fmt::Display> core::fmt::Display for DiscreteInverseCdfError<T> {
+    #[cfg_attr(coverage_nightly, coverage(off))]
+    fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result {
+        match self {
+            DiscreteInverseCdfError::ArgumentOutOfRange => {
+                write!(f, "argument is outside [0, 1]")
+            }
+            DiscreteInverseCdfError::NotRepresentable(v) => {
+                write!(f, "quantile ({v}) is not representable in the variate type")
+            }
+        }
+    }
+}
+
+impl<T: core::fmt::Debug + core::fmt::Display> core::error::Error for DiscreteInverseCdfError<T> {}
+
 /// The `ContinuousCDF` trait is used to specify an interface for univariate
 /// distributions for which cdf float arguments are sensible.
 pub trait ContinuousCDF<K: Float, T: Float>: Min<K> + Max<K> {
@@ -224,7 +257,11 @@ pub trait ContinuousCDF<K: Float, T: Float>: Min<K> + Max<K> {
     #[doc(alias = "quantile function")]
     #[doc(alias = "quantile")]
     fn try_inverse_cdf(&self, p: T) -> Result<K, InverseCdfError> {
-        Ok(self.inverse_cdf(p))
+        if !(T::zero()..=T::one()).contains(&p) {
+            Err(InverseCdfError::ArgumentOutOfRange)
+        } else {
+            Ok(self.inverse_cdf(p))
+        }
     }
 }
 
@@ -265,15 +302,13 @@ pub trait DiscreteCDF<K: Sized + Num + Ord + Clone + NumAssignOps, T: Float>:
     /// Due to issues with rounding and floating-point accuracy the default implementation may be ill-behaved
     /// Specialized inverse cdfs should be used whenever possible.
     ///
-    /// # Panics
-    /// this default impl panics if provided `p` not on interval [0.0, 1.0]
+    /// Does not check that `p` lies on `[0.0, 1.0]`; use [`try_inverse_cdf`](DiscreteCDF::try_inverse_cdf)
+    /// if `p` is not already known to be valid.
     fn inverse_cdf(&self, p: T) -> K {
         if p <= self.cdf(self.min()) {
             return self.min();
         } else if p == T::one() {
             return self.max();
-        } else if !(T::zero()..=T::one()).contains(&p) {
-            panic!("p must be on [0, 1]")
         }
 
         let two = K::one() + K::one();
@@ -285,6 +320,17 @@ pub trait DiscreteCDF<K: Sized + Num + Ord + Clone + NumAssignOps, T: Float>:
 
         internal::integral_bisection_search(|p| self.cdf(p.clone()), p, lb, ub).unwrap()
     }
+
+    /// Due to issues with rounding and floating-point accuracy the default
+    /// implementation may be ill-behaved.
+    /// Specialized inverse cdfs should be used whenever possible.
+    fn try_inverse_cdf(&self, p: T) -> Result<K, DiscreteInverseCdfError<T>> {
+        if !(T::zero()..=T::one()).contains(&p) {
+            Err(DiscreteInverseCdfError::ArgumentOutOfRange)
+        } else {
+            Ok(self.inverse_cdf(p))
+        }
+    }
 }
 
 /// The `Continuous` trait  provides an interface for interacting with
diff --git a/src/distribution/normal.rs b/src/distribution/normal.rs
index d555cef8..bfc81940 100644
--- a/src/distribution/normal.rs
+++ b/src/distribution/normal.rs
@@ -159,10 +159,6 @@ impl ContinuousCDF<f64, f64> for Normal {
     /// normal distribution at `x`.
     /// In other languages, such as R, this is known as the the quantile function.
     ///
-    /// # Panics
-    ///
-    /// If `x < 0.0` or `x > 1.0`
-    ///
     /// # Formula
     ///
     /// ```text
@@ -172,11 +168,7 @@ impl ContinuousCDF<f64, f64> for Normal {
     /// where `μ` is the mean, `σ` is the standard deviation and `erfc_inv` is
     /// the inverse of the complementary error function
     fn inverse_cdf(&self, x: f64) -> f64 {
-        if !(0.0..=1.0).contains(&x) {
-            panic!("x must be in [0, 1]");
-        } else {
-            self.mean - (self.std_dev * f64_consts::SQRT_2 * erf::erfc_inv(2.0 * x))
-        }
+        self.mean - (self.std_dev * f64_consts::SQRT_2 * erf::erfc_inv(2.0 * x))
     }
 }
 
diff --git a/src/distribution/pareto.rs b/src/distribution/pareto.rs
index 241a240a..9f924c7f 100644
--- a/src/distribution/pareto.rs
+++ b/src/distribution/pareto.rs
@@ -181,11 +181,7 @@ impl ContinuousCDF<f64, f64> for Pareto {
     ///
     /// where `x_m` is the scale and `α` is the shape
     fn inverse_cdf(&self, p: f64) -> f64 {
-        if !(0.0..=1.0).contains(&p) {
-            panic!("x must be in [0, 1]");
-        } else {
-            self.scale * (1.0 - p).powf(-1.0 / self.shape)
-        }
+        self.scale * (1.0 - p).powf(-1.0 / self.shape)
     }
 }
 
diff --git a/src/distribution/students_t.rs b/src/distribution/students_t.rs
index d1130c65..1743f666 100644
--- a/src/distribution/students_t.rs
+++ b/src/distribution/students_t.rs
@@ -220,7 +220,6 @@ impl ContinuousCDF<f64, f64> for StudentsT {
     /// Student's T-distribution at `x`
     fn inverse_cdf(&self, x: f64) -> f64 {
         // first calculate inverse_cdf for normal Student's T
-        assert!((0.0..=1.0).contains(&x));
         let x1 = if x >= 0.5 { 1.0 - x } else { x };
         let a = 0.5 * self.freedom;
         let b = 0.5;
diff --git a/src/distribution/triangular.rs b/src/distribution/triangular.rs
index 8f926116..73ac85a9 100644
--- a/src/distribution/triangular.rs
+++ b/src/distribution/triangular.rs
@@ -240,10 +240,6 @@ impl ContinuousCDF<f64, f64> for Triangular {
         let a = self.min;
         let b = self.max;
         let c = self.mode;
-        if !(0.0..=1.0).contains(&p) {
-            panic!("x must be in [0, 1]");
-        }
-
         if p < (c - a) / (b - a) {
             a + ((c - a) * (b - a) * p).sqrt()
         } else {
diff --git a/src/distribution/uniform.rs b/src/distribution/uniform.rs
index 07b7c503..8a655a21 100644
--- a/src/distribution/uniform.rs
+++ b/src/distribution/uniform.rs
@@ -197,9 +197,7 @@ impl ContinuousCDF<f64, f64> for Uniform {
 
     /// Finds the value of `x` where `F(p) = x`
     fn inverse_cdf(&self, p: f64) -> f64 {
-        if !(0.0..=1.0).contains(&p) {
-            panic!("p must be in [0, 1], was {p}");
-        } else if p == 0.0 {
+        if p == 0.0 {
             self.min
         } else if p == 1.0 {
             self.max
diff --git a/src/distribution/weibull.rs b/src/distribution/weibull.rs
index f61b163f..f3f89bca 100644
--- a/src/distribution/weibull.rs
+++ b/src/distribution/weibull.rs
@@ -179,10 +179,6 @@ impl ContinuousCDF<f64, f64> for Weibull {
     ///
     /// where `k` is the shape and `λ` is the scale
     fn inverse_cdf(&self, p: f64) -> f64 {
-        if !(0.0..=1.0).contains(&p) {
-            panic!("x must be in [0, 1]");
-        }
-
         (-((-p).ln_1p() / self.scale_pow_shape_inv)).powf(1.0 / self.shape)
     }
 }