Correct lgamma workaround for very small arguments.

This allows simplification of the student's T quantile code. Added some more lgamma tests.
surajkulriya · Jul 25, 2020 · a316c72 · a316c72
1 parent 2153dfa
commit a316c72
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 2 deletions.
diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp
@@ -697,7 +697,7 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
       }
       if (overflow || !(boost::math::isfinite)(bet))
       {
-         xg = exp((boost::math::lgamma(a + 1, pol) + boost::math::lgamma(b + 1, pol) - boost::math::lgamma(a + b + 1, pol) + log(p) - log(b) + log(a + b)) / a);
+         xg = exp((boost::math::lgamma(a + 1, pol) + boost::math::lgamma(b, pol) - boost::math::lgamma(a + b, pol) + log(p)) / a);
       }
       else
          xg = pow(a * p * bet, 1/a);

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
@@ -240,7 +240,7 @@ T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
    {
       if (0 == z)
          return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
-      if (fabs(z) < 1 / tools::max_value<T>())
+      if (4 * fabs(z) < tools::epsilon<T>())
          result = -log(fabs(z));
       else
          result = log(fabs(1 / z - constants::euler<T>()));

diff --git a/test/test_gamma.hpp b/test/test_gamma.hpp
@@ -333,5 +333,17 @@ void test_spots(T, const char* name)
       BOOST_CHECK_CLOSE(::boost::math::lgamma(ldexp(static_cast<T>(11103367432951928LL), 62)), static_cast<T>(4.0411767712186990905102512019058204792570873633363159e36L), tolerance);
       BOOST_CHECK_CLOSE(::boost::math::lgamma(ldexp(static_cast<T>(11103367432951928LL), 326)), static_cast<T>(3.9754720509185529233002820161357111676582583112671658e116L), tolerance);
    }
+   //
+   // Super small values may cause spurious overflow:
+   //
+   if (std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::has_denorm)
+   {
+      T value = (std::numeric_limits<T>::min)();
+      while (value != 0)
+      {
+         BOOST_CHECK((boost::math::isfinite)(boost::math::lgamma(value)));
+         value /= 2;
+      }
+   }
 }