diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 458dc40a3d..e67752861d 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1290,7 +1290,6 @@ BOOST_MATH_GPU_ENABLED T incomplete_tgamma_large_x(const T& a, const T& x, const return result; } - // // Main incomplete gamma entry point, handles all four incomplete gamma's: // @@ -1813,6 +1812,48 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_imp(T a, T x, const Policy& pol) return log(gamma_q(a, x, pol)); } +// Calculate log of incomplete gamma function +template +BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol) +{ + using namespace boost::math; // temporary until we're in the right namespace + + BOOST_MATH_STD_USING_CORE + + // Check for invalid inputs (a < 0 or x < 0) + constexpr auto function = "boost::math::lgamma_p<%1%>(%1%, %1%)"; + if(a <= 0) + return policies::raise_domain_error(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol); + if(x < 0) + return policies::raise_domain_error(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol); + + // Taken from conditions on Line 1709. There are also + // conditions on Line 1368, but didn't implement that one here. + // The second condition ensures floats do not return -inf for small + // values of x. + if (((a > 4 * x) && (a >= max_factorial::value)) || ((T(-0.4) / log(x) < a) && (x < T(1.0)))){ + return log(detail::lower_gamma_series(a, x, pol)) - log(a) + a * log(x) - x - lgamma(a, pol); + } + + // + // Can't do better than taking the log of P, but... + // + // Figure out whether we need P or Q, since if we calculate P and it's too close to unity + // we will lose precision in the result, selection logic here is extracted from gamma_incomplete_imp_final: + // + bool need_p = false; + if ((x < 0.5) && (T(-0.4) / log(x) < a)) + need_p = true; + else if ((x < 1.1) && (x >= 0.5) && (x * 0.75f < a)) + need_p = true; + else if ((x < a) && (x >= 1.1)) + need_p = true; + + if (need_p) + return log(gamma_p(a, x, pol)); + return log1p(-gamma_q(a, x, pol), pol); +} + // // Ratios of two gamma functions: // @@ -2454,6 +2495,29 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t lgamma_q(T1 a, T2 z) { return lgamma_q(a, z, policies::policy<>()); } + +template +BOOST_MATH_GPU_ENABLED inline tools::promote_args_t lgamma_p(T1 a, T2 z, const Policy& /* pol */) +{ + typedef tools::promote_args_t result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + return policies::checked_narrowing_cast( + detail::lgamma_incomplete_lower_imp(static_cast(a), + static_cast(z), forwarding_policy()), "lgamma_p<%1%>(%1%, %1%)"); +} + +template +BOOST_MATH_GPU_ENABLED inline tools::promote_args_t lgamma_p(T1 a, T2 z) +{ + return lgamma_p(a, z, policies::policy<>()); +} // // Regularised lower incomplete gamma: // diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp index ca16fc8d96..6c8394ddf9 100644 --- a/include/boost/math/special_functions/math_fwd.hpp +++ b/include/boost/math/special_functions/math_fwd.hpp @@ -567,6 +567,12 @@ namespace boost template BOOST_MATH_GPU_ENABLED tools::promote_args_t lgamma_q(RT1 a, RT2 z, const Policy&); + template + BOOST_MATH_GPU_ENABLED tools::promote_args_t lgamma_p(RT1 a, RT2 z); + + template + BOOST_MATH_GPU_ENABLED tools::promote_args_t lgamma_p(RT1 a, RT2 z, const Policy&); + template BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_p(RT1 a, RT2 z); @@ -1525,6 +1531,9 @@ namespace boost \ template \ BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t lgamma_q(RT1 a, RT2 z){ return boost::math::lgamma_q(a, z, Policy()); }\ +\ + template \ + BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t lgamma_p(RT1 a, RT2 z){ return boost::math::lgamma_p(a, z, Policy()); }\ \ template \ BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t gamma_p(RT1 a, RT2 z){ return boost::math::gamma_p(a, z, Policy()); }\ diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp index 88f0bca174..b12fe8568c 100644 --- a/test/test_igamma.hpp +++ b/test/test_igamma.hpp @@ -263,7 +263,18 @@ void test_spots(T, const char* name = nullptr) BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast(501.25), static_cast(2000)), static_cast(-810.2453406781655559126505101822969531699112391075198076300675402L), tolerance); BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast(20), static_cast(0.25)), static_cast(-2.946458104491857816330873290969917497748067639461638294404e-31L), tolerance); BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast(40), static_cast(0.75)), static_cast(-5.930604927955460343652485525435087275997461623988991819824e-54L), tolerance); -#if defined(__CYGWIN__) || defined(__MINGW32__) + + // + // Check that lgamma_q returns correct values with spot values calculated via wolframalpha log(P[a, x]) + // + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(500), static_cast(10)), static_cast(-1470.017750815998931281954666549641187420649099004671023115157832L), tolerance); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(100), static_cast(0.10)), static_cast(-594.0968942751157965882143339808498054972378571169718495813176479L), 2*tolerance); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(20), static_cast(10.25)), static_cast(-5.404004887981642339930593767572610169901594898478031307722239712L), tolerance); + // Small "a" produce larger errors + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(0.1), static_cast(100)), static_cast(-6.142676307812300286394850046780157761360272418143933840132692383e-47L), 20*tolerance); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(0.025), static_cast(10)), static_cast(-0.0000001118627574013195004426530864010428145980391578084985091741192526L), 15*tolerance); + + #if defined(__CYGWIN__) || defined(__MINGW32__) T gcc_win_mul = 2; #else T gcc_win_mul = 1; @@ -287,6 +298,26 @@ void test_spots(T, const char* name = nullptr) BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast(1200), static_cast(1250.25)), static_cast(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits::max_digits10 >= 36 || std::is_same::value) ? 750 : (std::is_same::value ? 1 : 50))); // Test fails on ARM64 and s390x long doubles and real_concept types unless tolerance is adjusted BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast(2200), static_cast(2249.75)), static_cast(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point::value ? 1 : 10)); BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast(2200), static_cast(2250.25)), static_cast(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same::value ? 1 : (std::is_floating_point::value ? 100 : 200))); + + // Long double and real_concept types need increased precision + T real_concept_tol = 1; + if (std::is_same::value || std::is_same::value){ + real_concept_tol = 3; + } + // Pair of tests that bisect the crossover condition in our code at double and then quad precision + // Oddly, the crossover condition is smaller for quad precision. This is because max_factorial is 100 + // for boost::multiprecision::cpp_bin_float_quad and 170 for doubles. + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(169.75), static_cast(0.6)), static_cast(-792.5976910460276693711943228399034564720784916702808816810639048L), tolerance * real_concept_tol); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(170.25), static_cast(0.6)), static_cast(-795.4224823388425178061598242301088759551669254174241556830902684L), tolerance * real_concept_tol); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(99.75), static_cast(0.6)), static_cast(-414.1360294778432398276991161095707673136841267685069121354421888L), tolerance * real_concept_tol); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(100.25), static_cast(0.6)), static_cast(-416.6965522800555653440139072333570825784615163822944047186474994L), tolerance * real_concept_tol); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(43.7796314883057249289777246303856372833251953125), static_cast(0.99)), static_cast(-125.889061559356008501708856783807277679443359375L), tolerance * real_concept_tol); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(39.799664989368835676941671408712863922119140625), static_cast(0.99)), static_cast(-110.9450811443740741159224346421309229458874083101891810411634707L), tolerance * real_concept_tol); + + // Check large a, x values. Precision just isn't great here. + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(1450.25), static_cast(1500.75)), static_cast(-0.09812447528127799786140178403478691390753413399549580160096975713L), tolerance * (std::is_same::value ? 8 : 1)); + BOOST_CHECK_CLOSE(::boost::math::lgamma_p(static_cast(2000), static_cast(1900)), static_cast(-4.448523733381445722945397105917814000790587922314824687065050805L), tolerance * gcc_win_mul * (std::is_same::value ? 8 : 1)); + // // Coverage: // @@ -302,6 +333,10 @@ void test_spots(T, const char* name = nullptr) BOOST_CHECK_THROW(boost::math::lgamma_q(static_cast(1), static_cast(-2)), std::domain_error); BOOST_CHECK_THROW(boost::math::lgamma_q(static_cast(0), static_cast(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::lgamma_p(static_cast(-1), static_cast(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::lgamma_p(static_cast(1), static_cast(-2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::lgamma_p(static_cast(0), static_cast(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast(-1), static_cast(2)), std::domain_error); BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast(1), static_cast(-2)), std::domain_error); BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast(0), static_cast(2)), std::domain_error); @@ -317,6 +352,10 @@ void test_spots(T, const char* name = nullptr) BOOST_CHECK((boost::math::isnan)(boost::math::lgamma_q(static_cast(1), static_cast(-2)))); BOOST_CHECK((boost::math::isnan)(boost::math::lgamma_q(static_cast(0), static_cast(2)))); + BOOST_CHECK((boost::math::isnan)(boost::math::lgamma_p(static_cast(-1), static_cast(2)))); + BOOST_CHECK((boost::math::isnan)(boost::math::lgamma_p(static_cast(1), static_cast(-2)))); + BOOST_CHECK((boost::math::isnan)(boost::math::lgamma_p(static_cast(0), static_cast(2)))); + BOOST_CHECK((boost::math::isnan)(boost::math::gamma_p_derivative(static_cast(-1), static_cast(2)))); BOOST_CHECK((boost::math::isnan)(boost::math::gamma_p_derivative(static_cast(1), static_cast(-2)))); BOOST_CHECK((boost::math::isnan)(boost::math::gamma_p_derivative(static_cast(0), static_cast(2))));