diff --git a/include/boost/math/special_functions/detail/lanczos_sse2.hpp b/include/boost/math/special_functions/detail/lanczos_sse2.hpp index df1a04743..2ebde18e0 100644 --- a/include/boost/math/special_functions/detail/lanczos_sse2.hpp +++ b/include/boost/math/special_functions/detail/lanczos_sse2.hpp @@ -51,6 +51,15 @@ inline double lanczos13m53::lanczos_sum(const double& x) static_cast(23531376880.41075968857200767445163675473L), static_cast(0u) }; + + static const double lim = 4.31965e+25; // By experiment, the largest x for which the SSE2 code does not go bad. + + if (x > lim) + { + double z = 1 / x; + return ((((((((((((coeff[24] * z + coeff[22]) * z + coeff[20]) * z + coeff[18]) * z + coeff[16]) * z + coeff[14]) * z + coeff[12]) * z + coeff[10]) * z + coeff[8]) * z + coeff[6]) * z + coeff[4]) * z + coeff[2]) * z + coeff[0]) / ((((((((((((coeff[25] * z + coeff[23]) * z + coeff[21]) * z + coeff[19]) * z + coeff[17]) * z + coeff[15]) * z + coeff[13]) * z + coeff[11]) * z + coeff[9]) * z + coeff[7]) * z + coeff[5]) * z + coeff[3]) * z + coeff[1]); + } + __m128d vx = _mm_load1_pd(&x); __m128d sum_even = _mm_load_pd(coeff); __m128d sum_odd = _mm_load_pd(coeff+2); @@ -136,6 +145,15 @@ inline double lanczos13m53::lanczos_sum_expG_scaled(const double& x) static_cast(56906521.91347156388090791033559122686859L), static_cast(0u) }; + + static const double lim = 4.76886e+25; // By experiment, the largest x for which the SSE2 code does not go bad. + + if (x > lim) + { + double z = 1 / x; + return ((((((((((((coeff[24] * z + coeff[22]) * z + coeff[20]) * z + coeff[18]) * z + coeff[16]) * z + coeff[14]) * z + coeff[12]) * z + coeff[10]) * z + coeff[8]) * z + coeff[6]) * z + coeff[4]) * z + coeff[2]) * z + coeff[0]) / ((((((((((((coeff[25] * z + coeff[23]) * z + coeff[21]) * z + coeff[19]) * z + coeff[17]) * z + coeff[15]) * z + coeff[13]) * z + coeff[11]) * z + coeff[9]) * z + coeff[7]) * z + coeff[5]) * z + coeff[3]) * z + coeff[1]); + } + __m128d vx = _mm_load1_pd(&x); __m128d sum_even = _mm_load_pd(coeff); __m128d sum_odd = _mm_load_pd(coeff+2); diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 417a67426..04dd71316 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -277,7 +277,11 @@ T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0) T zgh = static_cast(z + Lanczos::g() - boost::math::constants::half()); result = log(zgh) - 1; result *= z - 0.5f; - result += log(Lanczos::lanczos_sum_expG_scaled(z)); + // + // Only add on the lanczos sum part if we're going to need it: + // + if(result * tools::epsilon() < 20) + result += log(Lanczos::lanczos_sum_expG_scaled(z)); } if(sign) diff --git a/test/test_gamma.hpp b/test/test_gamma.hpp index f259e963a..47cd3d5eb 100644 --- a/test/test_gamma.hpp +++ b/test/test_gamma.hpp @@ -324,5 +324,14 @@ void test_spots(T, const char* name) BOOST_CHECK_EQUAL(boost::math::tgamma(-std::numeric_limits::denorm_min()), -std::numeric_limits::infinity()); BOOST_CHECK_EQUAL(boost::math::tgamma(std::numeric_limits::denorm_min()), std::numeric_limits::infinity()); } + // + // Extra large values for lgamma, see https://github.com/boostorg/math/issues/242 + // + if (boost::math::tools::digits() >= std::numeric_limits::digits) + { + BOOST_CHECK_CLOSE(::boost::math::lgamma(ldexp(static_cast(11103367432951928LL), 32)), static_cast(2.7719825960021351251696385101478518546793793286704974382373670822285114741208958e27L), tolerance); + BOOST_CHECK_CLOSE(::boost::math::lgamma(ldexp(static_cast(11103367432951928LL), 62)), static_cast(4.0411767712186990905102512019058204792570873633363159e36L), tolerance); + BOOST_CHECK_CLOSE(::boost::math::lgamma(ldexp(static_cast(11103367432951928LL), 326)), static_cast(3.9754720509185529233002820161357111676582583112671658e116), tolerance); + } }