fix: Finish large beta edge handling

Author: day01

src/function/beta.rs

@@ -3,11 +3,10 @@
 //!
 //! This module sets the default precision more tightly than crate defaults for `DEFAULT_EPS`
 
-mod double_double;
 mod large_params;
 mod temme;
 
-use crate::function::gamma;
+use crate::function::{double_double, gamma};
 use crate::prec;
 #[cfg(not(feature = "std"))]
 use num_traits::Float as _;
@@ -295,6 +294,14 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result<f64, BetaFuncError> {
         return Ok(x);
     }
 
+    if b == 1.0 && a + b == a {
+        return Ok((a * x.ln()).exp());
+    }
+
+    if a == 1.0 && a + b == b {
+        return Ok(-(b * (-x).ln_1p()).exp_m1());
+    }
+
     let log_prefactor = match large_params::log_prefactor(a, b, x) {
         Some(large_params::LogPrefactor::Value(logarithm)) => logarithm,
         Some(large_params::LogPrefactor::Underflow) => {
@@ -706,6 +713,16 @@ mod tests {
         assert_eq!(beta_reg(1e8, 2e8, 1.0), 1.0);
     }
 
+    #[test]
+    fn test_beta_reg_extreme_shape_one() {
+        assert_eq!(checked_beta_reg(1e308, 1.0, 0.0), Ok(0.0));
+        assert_eq!(checked_beta_reg(1e308, 1.0, 0.5), Ok(0.0));
+        assert_eq!(checked_beta_reg(1e308, 1.0, 1.0), Ok(1.0));
+        assert_eq!(checked_beta_reg(1.0, 1e308, 0.0), Ok(0.0));
+        assert_eq!(checked_beta_reg(1.0, 1e308, 0.5), Ok(1.0));
+        assert_eq!(checked_beta_reg(1.0, 1e308, 1.0), Ok(1.0));
+    }
+
     #[test]
     fn test_beta_reg_next_down_from_one_is_not_a_boundary() {
         let x = f64::from_bits(1.0_f64.to_bits() - 1);

src/function/beta/large_params.rs

@@ -1,5 +1,5 @@
-use super::double_double::{accurate_ln_one_plus, add, divide, multiply, two_sum};
 use crate::consts;
+use crate::function::double_double::{accurate_ln_one_plus, add, divide, multiply, two_sum};
 #[cfg(not(feature = "std"))]
 use num_traits::Float as _;
 

src/function/beta/temme.rs

@@ -1,6 +1,6 @@
-use super::double_double::{add, divide, exp, multiply, two_sum};
 use super::large_params::log_ratio;
 use crate::consts;
+use crate::function::double_double::{add, divide, exp, multiply, two_sum};
 use crate::function::erf;
 #[cfg(not(feature = "std"))]
 use num_traits::Float as _;

src/function/beta/double_double.rs src/function/double_double.rssrc/function/beta/double_double.rs renamed to src/function/double_double.rs

@@ -1,6 +1,8 @@
 //! Provides a host of special statistical functions (e.g. the beta function or
 //! the error function)
 
+mod double_double;
+
 pub mod beta;
 pub mod erf;
 pub mod evaluate;