Avoid panic on large alpha for NormalInverseGaussian

The new calculation method avoids calculating alpha*alpha, which could
produce `inf` and lead to invalid parameters for InverseGaussian.

Also, an explicit error is now returned when alpha is inf.

The new parameter calculation can slightly affect the values sampled
by the distribution.

Author: mstoeckl

CHANGELOG.md

@@ -19,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixes
 - Fix `Geometric::new` for small `p > 0` where `1 - p` rounds to 1 (#36)
 - Fix panic in `FisherF::new` on almost zero parameters (#39)
+- Fix panic in `NormalInverseGaussian::new` with very large `alpha`; this is a Value-breaking change (#40)
 
 ## [0.5.2]
 

src/normal_inverse_gaussian.rs

@@ -1,4 +1,4 @@
-use crate::{Distribution, InverseGaussian, StandardNormal, StandardUniform};
+use crate::{Distribution, InverseGaussian, InverseGaussianError, StandardNormal, StandardUniform};
 use core::fmt;
 use num_traits::Float;
 use rand::Rng;
@@ -8,6 +8,8 @@ use rand::Rng;
 pub enum Error {
     /// `alpha <= 0` or `nan`.
     AlphaNegativeOrNull,
+    /// `alpha` is `inf` (or, if subnormals are disabled, too close to the maximum finite float value)
+    AlphaInfinite,
     /// `|beta| >= alpha` or `nan`.
     AbsoluteBetaNotLessThanAlpha,
 }
@@ -18,6 +20,9 @@ impl fmt::Display for Error {
             Error::AlphaNegativeOrNull => {
                 "alpha <= 0 or is NaN in normal inverse Gaussian distribution"
             }
+            Error::AlphaInfinite => {
+                "alpha is +infinity (or too close to the maximum finite value, if subnormal numbers are not supported) in normal inverse Gaussian distribution"
+            }
             Error::AbsoluteBetaNotLessThanAlpha => {
                 "|beta| >= alpha or is NaN in normal inverse Gaussian distribution"
             }
@@ -68,6 +73,9 @@ where
 {
     /// Construct a new `NormalInverseGaussian` distribution with the given alpha (tail heaviness) and
     /// beta (asymmetry) parameters.
+    ///
+    /// Note: If subnormal numbers are not supported or are disabled, this sampler may panic or produce
+    /// incorrect output if alpha is close to the maximum finite float value.
     pub fn new(alpha: F, beta: F) -> Result<NormalInverseGaussian<F>, Error> {
         if !(alpha > F::zero()) {
             return Err(Error::AlphaNegativeOrNull);
@@ -76,12 +84,16 @@ where
         if !(beta.abs() < alpha) {
             return Err(Error::AbsoluteBetaNotLessThanAlpha);
         }
-
-        let gamma = (alpha * alpha - beta * beta).sqrt();
-
+        // Note: this calculation method for gamma = sqrt(alpha * alpha - beta * beta)
+        // avoids overflow if alpha is large, ensuring gamma <= alpha, which implies
+        // (assuming IEEE754 with subnormals) mu = 1.0 / gamma >= 1 / F::max_value() > 0.
+        let r = beta / alpha;
+        let gamma = alpha * (F::one() - r * r).sqrt();
         let mu = F::one() / gamma;
-
-        let inverse_gaussian = InverseGaussian::new(mu, F::one()).unwrap();
+        let inverse_gaussian = InverseGaussian::new(mu, F::one()).map_err(|x| match x {
+            InverseGaussianError::MeanNegativeOrNull => Error::AlphaInfinite,
+            InverseGaussianError::ShapeNegativeOrNull => unreachable!(),
+        })?;
 
         Ok(Self {
             beta,
@@ -124,6 +136,7 @@ mod tests {
         assert!(NormalInverseGaussian::new(-1.0, 1.0).is_err());
         assert!(NormalInverseGaussian::new(-1.0, -1.0).is_err());
         assert!(NormalInverseGaussian::new(1.0, 2.0).is_err());
+        assert!(NormalInverseGaussian::new(f64::INFINITY, 2.0).is_err());
         assert!(NormalInverseGaussian::new(2.0, 1.0).is_ok());
     }