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.

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` (#40)
 
 ## [0.5.2]
 

src/normal_inverse_gaussian.rs

@@ -8,6 +8,8 @@ use rand::Rng;
 pub enum Error {
     /// `alpha <= 0` or `nan`.
     AlphaNegativeOrNull,
+    /// `alpha` is `inf`
+    AlphaInfinite,
     /// `|beta| >= alpha` or `nan`.
     AbsoluteBetaNotLessThanAlpha,
 }
@@ -18,6 +20,7 @@ impl fmt::Display for Error {
             Error::AlphaNegativeOrNull => {
                 "alpha <= 0 or is NaN in normal inverse Gaussian distribution"
             }
+            Error::AlphaInfinite => "alpha is +infinity in normal inverse Gaussian distribution",
             Error::AbsoluteBetaNotLessThanAlpha => {
                 "|beta| >= alpha or is NaN in normal inverse Gaussian distribution"
             }
@@ -73,14 +76,19 @@ where
             return Err(Error::AlphaNegativeOrNull);
         }
 
+        if alpha == F::infinity() {
+            return Err(Error::AlphaInfinite);
+        }
+
         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();
 
         Ok(Self {