2503 lines · plain
1// Copyright John Maddock 2006-7, 2013-20.2// Copyright Paul A. Bristow 2007, 2013-14.3// Copyright Nikhar Agrawal 2013-144// Copyright Christopher Kormanyos 2013-14, 2020, 20245// Copyright Matt Borland 2024.6// Use, modification and distribution are subject to the7// Boost Software License, Version 1.0. (See accompanying file8// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)9 10#ifndef BOOST_MATH_SF_GAMMA_HPP11#define BOOST_MATH_SF_GAMMA_HPP12 13#ifdef _MSC_VER14#pragma once15#endif16 17#include <boost/math/tools/config.hpp>18#include <boost/math/tools/series.hpp>19#include <boost/math/tools/fraction.hpp>20#include <boost/math/tools/precision.hpp>21#include <boost/math/tools/promotion.hpp>22#include <boost/math/tools/type_traits.hpp>23#include <boost/math/tools/numeric_limits.hpp>24#include <boost/math/tools/cstdint.hpp>25#include <boost/math/tools/assert.hpp>26#include <boost/math/policies/error_handling.hpp>27#include <boost/math/constants/constants.hpp>28#include <boost/math/special_functions/math_fwd.hpp>29#include <boost/math/special_functions/log1p.hpp>30#include <boost/math/special_functions/trunc.hpp>31#include <boost/math/special_functions/powm1.hpp>32#include <boost/math/special_functions/sqrt1pm1.hpp>33#include <boost/math/special_functions/lanczos.hpp>34#include <boost/math/special_functions/fpclassify.hpp>35#include <boost/math/special_functions/detail/igamma_large.hpp>36#include <boost/math/special_functions/detail/unchecked_factorial.hpp>37#include <boost/math/special_functions/detail/lgamma_small.hpp>38 39// Only needed for types larger than double40#ifndef BOOST_MATH_HAS_GPU_SUPPORT41#include <boost/math/special_functions/bernoulli.hpp>42#include <boost/math/special_functions/polygamma.hpp>43#endif44 45#ifdef _MSC_VER46# pragma warning(push)47# pragma warning(disable: 4702) // unreachable code (return after domain_error throw).48# pragma warning(disable: 4127) // conditional expression is constant.49# pragma warning(disable: 4100) // unreferenced formal parameter.50# pragma warning(disable: 6326) // potential comparison of a constant with another constant51// Several variables made comments,52// but some difficulty as whether referenced on not may depend on macro values.53// So to be safe, 4100 warnings suppressed.54// TODO - revisit this?55#endif56 57namespace boost{ namespace math{58 59namespace detail{60 61template <class T>62BOOST_MATH_GPU_ENABLED inline bool is_odd(T v, const boost::math::true_type&)63{64 int i = static_cast<int>(v);65 return i&1;66}67template <class T>68BOOST_MATH_GPU_ENABLED inline bool is_odd(T v, const boost::math::false_type&)69{70 // Oh dear can't cast T to int!71 BOOST_MATH_STD_USING72 T modulus = v - 2 * floor(v/2);73 return static_cast<bool>(modulus != 0);74}75template <class T>76BOOST_MATH_GPU_ENABLED inline bool is_odd(T v)77{78 return is_odd(v, ::boost::math::is_convertible<T, int>());79}80 81template <class T>82BOOST_MATH_GPU_ENABLED T sinpx(T z)83{84 // Ad hoc function calculates x * sin(pi * x),85 // taking extra care near when x is near a whole number.86 BOOST_MATH_STD_USING87 int sign = 1;88 if(z < 0)89 {90 z = -z;91 }92 T fl = floor(z);93 T dist; // LCOV_EXCL_LINE94 if(is_odd(fl))95 {96 fl += 1;97 dist = fl - z;98 sign = -sign;99 }100 else101 {102 dist = z - fl;103 }104 BOOST_MATH_ASSERT(fl >= 0);105 if(dist > T(0.5))106 dist = 1 - dist;107 T result = sin(dist*boost::math::constants::pi<T>());108 return sign*z*result;109} // template <class T> T sinpx(T z)110//111// tgamma(z), with Lanczos support:112//113template <class T, class Policy, class Lanczos>114BOOST_MATH_GPU_ENABLED T gamma_imp_final(T z, const Policy& pol, const Lanczos& l)115{116 BOOST_MATH_STD_USING117 118 (void)l; // Suppresses unused variable warning when BOOST_MATH_INSTRUMENT is not defined119 120 T result = 1; 121 122#ifdef BOOST_MATH_INSTRUMENT123 static bool b = false;124 if(!b)125 {126 std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;127 b = true;128 }129#endif130 constexpr auto function = "boost::math::tgamma<%1%>(%1%)";131 132 if(z <= 0)133 {134 // shift z to > 1:135 while(z < 0)136 {137 result /= z;138 z += 1;139 }140 }141 BOOST_MATH_INSTRUMENT_VARIABLE(result);142 if((floor(z) == z) && (z < max_factorial<T>::value))143 {144 result *= unchecked_factorial<T>(static_cast<unsigned>(itrunc(z, pol) - 1));145 BOOST_MATH_INSTRUMENT_VARIABLE(result);146 }147 else if (z < tools::root_epsilon<T>())148 {149 if (z < 1 / tools::max_value<T>())150 result = policies::raise_overflow_error<T>(function, nullptr, pol);151 result *= 1 / z - constants::euler<T>();152 }153 else154 {155 result *= Lanczos::lanczos_sum(z);156 T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());157 T lzgh = log(zgh);158 BOOST_MATH_INSTRUMENT_VARIABLE(result);159 BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());160 if(z * lzgh > tools::log_max_value<T>())161 {162 // we're going to overflow unless this is done with care:163 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);164 if(lzgh * z / 2 > tools::log_max_value<T>())165 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);166 T hp = pow(zgh, T((z / 2) - T(0.25)));167 BOOST_MATH_INSTRUMENT_VARIABLE(hp);168 result *= hp / exp(zgh);169 BOOST_MATH_INSTRUMENT_VARIABLE(result);170 if(tools::max_value<T>() / hp < result)171 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);172 result *= hp;173 BOOST_MATH_INSTRUMENT_VARIABLE(result);174 }175 else176 {177 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);178 BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, T(z - boost::math::constants::half<T>())));179 BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));180 result *= pow(zgh, T(z - boost::math::constants::half<T>())) / exp(zgh);181 BOOST_MATH_INSTRUMENT_VARIABLE(result);182 }183 }184 return result;185}186 187#ifdef BOOST_MATH_ENABLE_CUDA188# pragma nv_diag_suppress 2190189#endif190 191// SYCL compilers can not support recursion so we extract it into a dispatch function192template <class T, class Policy, class Lanczos>193BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T gamma_imp(T z, const Policy& pol, const Lanczos& l)194{195 BOOST_MATH_STD_USING196 197 T result = 1;198 constexpr auto function = "boost::math::tgamma<%1%>(%1%)";199 200 if(z <= 0)201 {202 if(floor(z) == z)203 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);204 if(z <= -20)205 {206#ifndef BOOST_MATH_NO_EXCEPTIONS207 try208#endif209 {210 result = gamma_imp_final(T(-z), pol, l) * sinpx(z);211 }212#ifndef BOOST_MATH_NO_EXCEPTIONS213 catch (const std::overflow_error&)214 {215 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);216 }217#endif218 BOOST_MATH_INSTRUMENT_VARIABLE(result);219 BOOST_MATH_IF_CONSTEXPR(!boost::math::numeric_limits<T>::is_specialized || (boost::math::numeric_limits<T>::digits > 64))220 {221 if ((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))222 {223 return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE MP only.224 225 }226 }227 else228 {229 // Result can never be small: tgamma[-z] is always larger than sinpx[z] is small.230 // Specifically, sinpx can never be larger than 1 / epsilon which is too small to231 // ever generate a value less than one for `result`, unless T has a truely232 // exceptional number of digits precision.233 BOOST_MATH_ASSERT((fabs(result) > 1) || (tools::max_value<T>() * fabs(result) > boost::math::constants::pi<T>()));234 }235 result = -boost::math::constants::pi<T>() / result;236 if (result == 0)237 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);238 /*239 * Result can never be subnormal as we have a value > 1 in the numerator:240 if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)241 return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);242 */243 BOOST_MATH_INSTRUMENT_VARIABLE(result);244 return result;245 }246 }247 248 return gamma_imp_final(T(z), pol, l);249}250 251#ifdef BOOST_MATH_ENABLE_CUDA252# pragma nv_diag_default 2190253#endif254 255//256// lgamma(z) with Lanczos support:257//258template <class T, class Policy, class Lanczos>259BOOST_MATH_GPU_ENABLED T lgamma_imp_final(T z, const Policy& pol, const Lanczos& l, int* sign = nullptr)260{261#ifdef BOOST_MATH_INSTRUMENT262 static bool b = false;263 if(!b)264 {265 std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;266 b = true;267 }268#endif269 270 BOOST_MATH_STD_USING271 272 constexpr auto function = "boost::math::lgamma<%1%>(%1%)";273 274 T result = 0;275 int sresult = 1;276 277 if (z < tools::root_epsilon<T>())278 {279 if (0 == z)280 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);281 if (4 * fabs(z) < tools::epsilon<T>())282 result = -log(fabs(z));283 else284 result = log(fabs(1 / z - constants::euler<T>()));285 if (z < 0)286 sresult = -1;287 }288 else if(z < 15)289 {290 typedef typename policies::precision<T, Policy>::type precision_type;291 typedef boost::math::integral_constant<int,292 precision_type::value <= 0 ? 0 :293 precision_type::value <= 64 ? 64 :294 precision_type::value <= 113 ? 113 : 0295 > tag_type;296 297 result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);298 }299 else if((z >= 3) && (z < 100) && (boost::math::numeric_limits<T>::max_exponent >= 1024))300 {301 // taking the log of tgamma reduces the error, no danger of overflow here:302 result = log(gamma_imp(z, pol, l));303 }304 else305 {306 // regular evaluation:307 T zgh = static_cast<T>(z + T(Lanczos::g()) - boost::math::constants::half<T>());308 result = log(zgh) - 1;309 result *= z - 0.5f;310 //311 // Only add on the lanczos sum part if we're going to need it:312 //313 if(result * tools::epsilon<T>() < 20)314 result += log(Lanczos::lanczos_sum_expG_scaled(z));315 }316 317 if(sign)318 *sign = sresult;319 return result;320}321 322#ifdef BOOST_MATH_ENABLE_CUDA323# pragma nv_diag_suppress 2190324#endif325 326template <class T, class Policy, class Lanczos>327BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = nullptr)328{329 BOOST_MATH_STD_USING330 331 if(z <= -tools::root_epsilon<T>())332 {333 constexpr auto function = "boost::math::lgamma<%1%>(%1%)";334 335 T result = 0;336 int sresult = 1;337 338 // reflection formula:339 if(floor(z) == z)340 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);341 342 T t = sinpx(z);343 z = -z;344 if(t < 0)345 {346 t = -t;347 }348 else349 {350 sresult = -sresult;351 }352 result = log(boost::math::constants::pi<T>()) - lgamma_imp_final(T(z), pol, l) - log(t);353 354 if(sign)355 {356 *sign = sresult;357 }358 359 return result;360 }361 else362 {363 return lgamma_imp_final(T(z), pol, l, sign);364 }365}366 367#ifdef BOOST_MATH_ENABLE_CUDA368# pragma nv_diag_default 2190369#endif370 371//372// Incomplete gamma functions follow:373//374template <class T>375struct upper_incomplete_gamma_fract376{377private:378 T z, a;379 int k;380public:381 typedef boost::math::pair<T,T> result_type;382 383 BOOST_MATH_GPU_ENABLED upper_incomplete_gamma_fract(T a1, T z1)384 : z(z1-a1+1), a(a1), k(0)385 {386 }387 388 BOOST_MATH_GPU_ENABLED result_type operator()()389 {390 ++k;391 z += 2;392 return result_type(k * (a - k), z);393 }394};395 396template <class T>397BOOST_MATH_GPU_ENABLED inline T upper_gamma_fraction(T a, T z, T eps)398{399 // Multiply result by z^a * e^-z to get the full400 // upper incomplete integral. Divide by tgamma(z)401 // to normalise.402 upper_incomplete_gamma_fract<T> f(a, z);403 return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));404}405 406template <class T>407struct lower_incomplete_gamma_series408{409private:410 T a, z, result;411public:412 typedef T result_type;413 BOOST_MATH_GPU_ENABLED lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}414 415 BOOST_MATH_GPU_ENABLED T operator()()416 {417 T r = result;418 a += 1;419 result *= z/a;420 return r;421 }422};423 424template <class T, class Policy>425BOOST_MATH_GPU_ENABLED inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)426{427 // Multiply result by ((z^a) * (e^-z) / a) to get the full428 // lower incomplete integral. Then divide by tgamma(a)429 // to get the normalised value.430 lower_incomplete_gamma_series<T> s(a, z);431 boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();432 T factor = policies::get_epsilon<T, Policy>();433 T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);434 policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);435 return result;436}437 438#ifndef BOOST_MATH_HAS_GPU_SUPPORT439 440//441// Fully generic tgamma and lgamma use Stirling's approximation442// with Bernoulli numbers.443//444template<class T>445boost::math::size_t highest_bernoulli_index()446{447 const float digits10_of_type = (boost::math::numeric_limits<T>::is_specialized448 ? static_cast<float>(boost::math::numeric_limits<T>::digits10)449 : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));450 451 // Find the high index n for Bn to produce the desired precision in Stirling's calculation.452 return static_cast<boost::math::size_t>(18.0F + (0.6F * digits10_of_type));453}454 455template<class T>456int minimum_argument_for_bernoulli_recursion()457{458 BOOST_MATH_STD_USING459 460 const float digits10_of_type = (boost::math::numeric_limits<T>::is_specialized461 ? (float) boost::math::numeric_limits<T>::digits10462 : (float) (boost::math::tools::digits<T>() * 0.301F));463 464 int min_arg = (int) (digits10_of_type * 1.7F);465 466 if(digits10_of_type < 50.0F)467 {468 // The following code sequence has been modified469 // within the context of issue 396.470 471 // The calculation of the test-variable limit has now472 // been protected against overflow/underflow dangers.473 474 // The previous line looked like this and did, in fact,475 // underflow ldexp when using certain multiprecision types.476 477 // const float limit = std::ceil(std::pow(1.0f / std::ldexp(1.0f, 1-boost::math::tools::digits<T>()), 1.0f / 20.0f));478 479 // The new safe version of the limit check is now here.480 const float d2_minus_one = ((digits10_of_type / 0.301F) - 1.0F);481 const float limit = ceil(exp((d2_minus_one * log(2.0F)) / 20.0F));482 483 min_arg = (int) (BOOST_MATH_GPU_SAFE_MIN(digits10_of_type * 1.7F, limit));484 }485 486 return min_arg;487}488 489template <class T, class Policy>490T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false)491{492 BOOST_MATH_STD_USING493 //494 // Calculates tgamma(z) / (z/e)^z495 // Requires that our argument is large enough for Sterling's approximation to hold.496 // Used internally when combining gamma's of similar magnitude without logarithms.497 //498 BOOST_MATH_ASSERT(minimum_argument_for_bernoulli_recursion<T>() <= z);499 500 // Perform the Bernoulli series expansion of Stirling's approximation.501 502 const boost::math::size_t number_of_bernoullis_b2n = policies::get_max_series_iterations<Policy>();503 504 T one_over_x_pow_two_n_minus_one = 1 / z;505 const T one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;506 T sum = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one;507 const T target_epsilon_to_break_loop = sum * boost::math::tools::epsilon<T>();508 const T half_ln_two_pi_over_z = sqrt(boost::math::constants::two_pi<T>() / z);509 T last_term = 2 * sum;510 511 for (boost::math::size_t n = 2U;; ++n)512 {513 one_over_x_pow_two_n_minus_one *= one_over_x2;514 515 const boost::math::size_t n2 = static_cast<boost::math::size_t>(n * 2U);516 517 const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U));518 519 if ((n >= 3U) && (abs(term) < target_epsilon_to_break_loop))520 {521 // We have reached the desired precision in Stirling's expansion.522 // Adding additional terms to the sum of this divergent asymptotic523 // expansion will not improve the result.524 525 // Break from the loop.526 break;527 }528 if (n > number_of_bernoullis_b2n)529 // Safety net, we hope to never get here:530 return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Exceeded maximum series iterations without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol); // LCOV_EXCL_LINE531 532 sum += term;533 534 // Sanity check for divergence:535 T fterm = fabs(term);536 if(fterm > last_term)537 // Safety net, we hope to never get here:538 return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Series became divergent without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol); // LCOV_EXCL_LINE539 last_term = fterm;540 }541 542 // Complete Stirling's approximation.543 T scaled_gamma_value = islog ? T(sum + log(half_ln_two_pi_over_z)) : T(exp(sum) * half_ln_two_pi_over_z);544 return scaled_gamma_value;545}546 547// Forward declaration of the lgamma_imp template specialization.548template <class T, class Policy>549T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = nullptr);550 551template <class T, class Policy>552T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)553{554 BOOST_MATH_STD_USING555 556 constexpr auto function = "boost::math::tgamma<%1%>(%1%)";557 558 // Check if the argument of tgamma is identically zero.559 const bool is_at_zero = (z == 0);560 561 if((boost::math::isnan)(z) || (is_at_zero) || ((boost::math::isinf)(z) && (z < 0)))562 return policies::raise_domain_error<T>(function, "Evaluation of tgamma at %1%.", z, pol);563 564 const bool b_neg = (z < 0);565 566 const bool floor_of_z_is_equal_to_z = (floor(z) == z);567 568 // Special case handling of small factorials:569 if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))570 {571 return boost::math::unchecked_factorial<T>(static_cast<unsigned>(itrunc(z) - 1));572 }573 574 // Make a local, unsigned copy of the input argument.575 T zz((!b_neg) ? z : -z);576 577 // Special case for ultra-small z:578 if(zz < tools::cbrt_epsilon<T>())579 {580 const T a0(1);581 const T a1(boost::math::constants::euler<T>());582 const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6);583 const T a2((six_euler_squared - boost::math::constants::pi_sqr<T>()) / 12);584 585 const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0);586 587 return 1 / inverse_tgamma_series;588 }589 590 // Scale the argument up for the calculation of lgamma,591 // and use downward recursion later for the final result.592 const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();593 594 int n_recur;595 596 if(zz < min_arg_for_recursion)597 {598 n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1;599 600 zz += n_recur;601 }602 else603 {604 n_recur = 0;605 }606 if (!n_recur)607 {608 if (zz > tools::log_max_value<T>())609 return b_neg ? policies::raise_underflow_error<T>(function, nullptr, pol) : policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE MP only610 if (log(zz) * zz / 2 > tools::log_max_value<T>())611 return b_neg ? policies::raise_underflow_error<T>(function, nullptr, pol) : policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE MP only612 }613 T gamma_value = scaled_tgamma_no_lanczos(zz, pol);614 T power_term = pow(zz, zz / 2);615 T exp_term = exp(-zz);616 gamma_value *= (power_term * exp_term);617 if (!n_recur && (tools::max_value<T>() / power_term < gamma_value))618 return b_neg ? policies::raise_underflow_error<T>(function, nullptr, pol) : policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE MP only619 gamma_value *= power_term;620 621 // Rescale the result using downward recursion if necessary.622 if(n_recur)623 {624 // The order of divides is important, if we keep subtracting 1 from zz625 // we DO NOT get back to z (cancellation error). Further if z < epsilon626 // we would end up dividing by zero. Also in order to prevent spurious627 // overflow with the first division, we must save dividing by |z| till last,628 // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z.629 zz = fabs(z) + 1;630 for(int k = 1; k < n_recur; ++k)631 {632 gamma_value /= zz;633 zz += 1;634 }635 gamma_value /= fabs(z);636 }637 638 // Return the result, accounting for possible negative arguments.639 if(b_neg)640 {641 // Provide special error analysis for:642 // * arguments in the neighborhood of a negative integer643 // * arguments exactly equal to a negative integer.644 645 // Check if the argument of tgamma is exactly equal to a negative integer.646 if(floor_of_z_is_equal_to_z)647 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); // LCOV_EXCL_LINE MP only648 649 T s = sinpx(z);650 if ((gamma_value > 1) && (tools::max_value<T>() / gamma_value < fabs(s)))651 return policies::raise_underflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE MP only652 gamma_value *= s; // LCOV_EXCL_LINE MP only653 654 BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);655 //656 // Result can never overflow, since sinpx(z) can never be smaller than machine epsilon and gamma_value > 1.657 //658 BOOST_MATH_ASSERT( (abs(gamma_value) > 1) || ((tools::max_value<T>() * abs(gamma_value)) > boost::math::constants::pi<T>())); // LCOV_EXCL_LINE MP only659 660 gamma_value = -boost::math::constants::pi<T>() / gamma_value;661 662 BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value); // LCOV_EXCL_LINE MP only663 //664 // We can never underflow since the numerator > 1 above and denominator is not infinite:665 //666 BOOST_MATH_ASSERT(gamma_value != 0); // LCOV_EXCL_LINE MP only667 }668 669 return gamma_value;670}671 672template <class T, class Policy>673inline T log_gamma_near_1(const T& z, Policy const& pol)674{675 //676 // This is for the multiprecision case where there is677 // no lanczos support, use a taylor series at z = 1,678 // see https://www.wolframalpha.com/input/?i=taylor+series+lgamma(x)+at+x+%3D+1679 //680 BOOST_MATH_STD_USING // ADL of std names681 682 // For some reason, several lines aren't triggered for coverage even though683 // adjacent lines are... weird!684 685 BOOST_MATH_ASSERT(fabs(z) < 1); // LCOV_EXCL_LINE686 687 T result = -constants::euler<T>() * z;688 689 T power_term = z * z / 2;690 int n = 2; // LCOV_EXCL_LINE691 T term = 0; // LCOV_EXCL_LINE692 693 do694 {695 term = power_term * boost::math::polygamma(n - 1, T(1), pol);696 result += term; // LCOV_EXCL_LINE697 ++n;698 power_term *= z / n;699 } while (fabs(result) * tools::epsilon<T>() < fabs(term));700 701 return result;702}703 704template <class T, class Policy>705T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)706{707 BOOST_MATH_STD_USING708 709 constexpr auto function = "boost::math::lgamma<%1%>(%1%)";710 711 // Check if the argument of lgamma is identically zero.712 const bool is_at_zero = (z == 0);713 714 if(is_at_zero)715 return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%.", z, pol);716 if((boost::math::isnan)(z))717 return policies::raise_domain_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);718 if((boost::math::isinf)(z))719 return policies::raise_overflow_error<T>(function, nullptr, pol);720 721 const bool b_neg = (z < 0);722 723 const bool floor_of_z_is_equal_to_z = (floor(z) == z);724 725 // Special case handling of small factorials:726 if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))727 {728 if (sign)729 *sign = 1; // LCOV_EXCL_LINE730 return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1));731 }732 733 // Make a local, unsigned copy of the input argument.734 T zz((!b_neg) ? z : -z);735 736 const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();737 738 T log_gamma_value;739 740 if (zz < min_arg_for_recursion)741 {742 // Here we simply take the logarithm of tgamma(). This is somewhat743 // inefficient, but simple. The rationale is that the argument here744 // is relatively small and overflow is not expected to be likely.745 if (sign)746 * sign = 1; // // LCOV_EXCL_LINE747 if(fabs(z - 1) < 0.25)748 {749 log_gamma_value = log_gamma_near_1(T(zz - 1), pol);750 }751 else if(fabs(z - 2) < 0.25)752 {753 log_gamma_value = log_gamma_near_1(T(zz - 2), pol) + log(zz - 1);754 }755 else if (z > -tools::root_epsilon<T>())756 {757 // Reflection formula may fail if z is very close to zero, let the series758 // expansion for tgamma close to zero do the work:759 if (sign)760 *sign = z < 0 ? -1 : 1; // LCOV_EXCL_LINE761 return log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos())));762 }763 else764 {765 // No issue with spurious overflow in reflection formula,766 // just fall through to regular code:767 T g = gamma_imp(zz, pol, lanczos::undefined_lanczos());768 if (sign)769 {770 *sign = g < 0 ? -1 : 1; // LCOV_EXCL_LINE MP only771 }772 log_gamma_value = log(abs(g));773 }774 }775 else776 {777 // Perform the Bernoulli series expansion of Stirling's approximation.778 T sum = scaled_tgamma_no_lanczos(zz, pol, true);779 log_gamma_value = zz * (log(zz) - 1) + sum;780 }781 782 int sign_of_result = 1;783 784 if(b_neg)785 {786 // Provide special error analysis if the argument is exactly787 // equal to a negative integer.788 789 // Check if the argument of lgamma is exactly equal to a negative integer.790 if(floor_of_z_is_equal_to_z)791 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol); // LCOV_EXCL_LINE MP only792 793 T t = sinpx(z);794 795 if(t < 0)796 {797 t = -t;798 }799 else800 {801 sign_of_result = -sign_of_result; // LCOV_EXCL_LINE MP only802 }803 804 log_gamma_value = - log_gamma_value + log(boost::math::constants::pi<T>()) - log(t);805 }806 807 if(sign != static_cast<int*>(nullptr)) { *sign = sign_of_result; }808 809 return log_gamma_value;810}811 812#endif // BOOST_MATH_HAS_GPU_SUPPORT813 814// In order for tgammap1m1_imp to compile we need a forward decl of boost::math::tgamma815// The rub is that we can't just use math_fwd so we provide one here only in that circumstance816#ifdef BOOST_MATH_HAS_NVRTC817template <class RT>818BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT> tgamma(RT z);819 820template <class RT1, class RT2>821BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> tgamma(RT1 a, RT2 z);822 823template <class RT1, class RT2, class Policy>824BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> tgamma(RT1 a, RT2 z, const Policy& pol);825#endif826 827//828// This helper calculates tgamma(dz+1)-1 without cancellation errors,829// used by the upper incomplete gamma with z < 1:830//831template <class T, class Policy, class Lanczos>832BOOST_MATH_GPU_ENABLED T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)833{834 BOOST_MATH_STD_USING835 836 typedef typename policies::precision<T,Policy>::type precision_type;837 838 typedef boost::math::integral_constant<int,839 precision_type::value <= 0 ? 0 :840 precision_type::value <= 64 ? 64 :841 precision_type::value <= 113 ? 113 : 0842 > tag_type;843 844 T result{};845 if(dz < 0)846 {847 if(dz < T(-0.5))848 {849 // Best method is simply to subtract 1 from tgamma:850 #ifdef BOOST_MATH_HAS_NVRTC851 result = ::tgamma(1+dz);852 #else853 result = boost::math::tgamma(1+dz, pol) - 1;854 #endif855 BOOST_MATH_INSTRUMENT_CODE(result);856 }857 else858 {859 // Use expm1 on lgamma:860 result = boost::math::expm1(-boost::math::log1p(dz, pol)861 + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l), pol);862 BOOST_MATH_INSTRUMENT_CODE(result);863 }864 }865 else866 {867 if(dz < 2)868 {869 // Use expm1 on lgamma:870 result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);871 BOOST_MATH_INSTRUMENT_CODE(result);872 }873 else874 {875 // Best method is simply to subtract 1 from tgamma:876 #ifdef BOOST_MATH_HAS_NVRTC877 result = ::tgamma(1+dz);878 #else879 result = boost::math::tgamma(1+dz, pol) - 1;880 #endif881 BOOST_MATH_INSTRUMENT_CODE(result);882 }883 }884 885 return result;886}887 888#ifndef BOOST_MATH_HAS_GPU_SUPPORT889 890template <class T, class Policy>891inline T tgammap1m1_imp(T z, Policy const& pol,892 const ::boost::math::lanczos::undefined_lanczos&)893{894 BOOST_MATH_STD_USING // ADL of std names895 896 if(fabs(z) < T(0.55))897 {898 return boost::math::expm1(log_gamma_near_1(z, pol));899 }900 return boost::math::expm1(boost::math::lgamma(1 + z, pol));901}902 903#endif // BOOST_MATH_HAS_GPU_SUPPORT904 905//906// Series representation for upper fraction when z is small:907//908template <class T>909struct small_gamma2_series910{911 typedef T result_type;912 913 BOOST_MATH_GPU_ENABLED small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}914 915 BOOST_MATH_GPU_ENABLED T operator()()916 {917 T r = result / (apn);918 result *= x;919 result /= ++n;920 apn += 1;921 return r;922 }923 924private:925 T result, x, apn;926 int n;927};928//929// calculate power term prefix (z^a)(e^-z) used in the non-normalised930// incomplete gammas:931//932template <class T, class Policy>933BOOST_MATH_GPU_ENABLED T full_igamma_prefix(T a, T z, const Policy& pol)934{935 BOOST_MATH_STD_USING936 937 if (z > tools::max_value<T>() || (a > 0 && z == 0))938 return 0;939 940 T alz = a * log(z);941 942 T prefix { };943 944 if(z >= 1)945 {946 if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))947 {948 prefix = pow(z, a) * exp(-z);949 }950 else if(a >= 1)951 {952 prefix = pow(T(z / exp(z/a)), a);953 }954 else955 {956 prefix = exp(alz - z); // LCOV_EXCL_LINE defensive programming, can probably never get here?957 }958 }959 else960 {961 if(alz > tools::log_min_value<T>())962 {963 prefix = pow(z, a) * exp(-z);964 }965 // LCOV_EXCL_START966 // Defensive programming, can probably never get here, very hard to prove though!967 else if(z/a < tools::log_max_value<T>())968 {969 prefix = pow(T(z / exp(z/a)), a);970 }971 else972 {973 prefix = exp(alz - z);974 }975 // LCOV_EXCL_STOP976 }977 //978 // This error handling isn't very good: it happens after the fact979 // rather than before it...980 // Typically though this method is used when the result is small, we should probably not overflow here...981 //982 if((boost::math::fpclassify)(prefix) == (int)BOOST_MATH_FP_INFINITE)983 return policies::raise_overflow_error<T>("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol); // LCOV_EXCL_LINE984 985 return prefix;986}987//988// Compute (z^a)(e^-z)/tgamma(a)989// most if the error occurs in this function:990//991template <class T, class Policy, class Lanczos>992BOOST_MATH_GPU_ENABLED T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)993{994 BOOST_MATH_STD_USING995 if (z >= tools::max_value<T>() || (a > 0 && z == 0))996 return 0;997 T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);998 T prefix{};999 T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;1000 1001 if(a < 1)1002 {1003 //1004 // We have to treat a < 1 as a special case because our Lanczos1005 // approximations are optimised against the factorials with a > 1,1006 // and for high precision types especially (128-bit reals for example)1007 // very small values of a can give rather erroneous results for gamma1008 // unless we do this:1009 //1010 // TODO: is this still required? Lanczos approx should be better now?1011 //1012 if((z <= tools::log_min_value<T>()) || (a < 1 / tools::max_value<T>()))1013 {1014 // Oh dear, have to use logs, should be free of cancellation errors though:1015 return exp(a * log(z) - z - lgamma_imp(a, pol, l));1016 }1017 else1018 {1019 // direct calculation, no danger of overflow as gamma(a) < 1/a1020 // for small a.1021 return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);1022 }1023 }1024 else if((fabs(d*d*a) <= 100) && (a > 150))1025 {1026 // special case for large a and a ~ z.1027 prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;1028 prefix = exp(prefix);1029 }1030 else1031 {1032 //1033 // general case.1034 // direct computation is most accurate, but use various fallbacks1035 // for different parts of the problem domain:1036 //1037 T alz = a * log(z / agh);1038 T amz = a - z;1039 if((BOOST_MATH_GPU_SAFE_MIN(alz, amz) <= tools::log_min_value<T>()) || (BOOST_MATH_GPU_SAFE_MAX(alz, amz) >= tools::log_max_value<T>()))1040 {1041 T amza = amz / a;1042 if((BOOST_MATH_GPU_SAFE_MIN(alz, amz)/2 > tools::log_min_value<T>()) && (BOOST_MATH_GPU_SAFE_MAX(alz, amz)/2 < tools::log_max_value<T>()))1043 {1044 // compute square root of the result and then square it:1045 T sq = pow(z / agh, a / 2) * exp(amz / 2);1046 prefix = sq * sq;1047 }1048 else if((BOOST_MATH_GPU_SAFE_MIN(alz, amz)/4 > tools::log_min_value<T>()) && (BOOST_MATH_GPU_SAFE_MAX(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))1049 {1050 // compute the 4th root of the result then square it twice:1051 T sq = pow(z / agh, a / 4) * exp(amz / 4);1052 prefix = sq * sq;1053 prefix *= prefix;1054 }1055 else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))1056 {1057 prefix = pow(T((z * exp(amza)) / agh), a);1058 }1059 else1060 {1061 prefix = exp(alz + amz);1062 }1063 }1064 else1065 {1066 prefix = pow(T(z / agh), a) * exp(amz);1067 }1068 }1069 prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);1070 return prefix;1071}1072 1073#ifndef BOOST_MATH_HAS_GPU_SUPPORT1074 1075//1076// And again, without Lanczos support:1077//1078template <class T, class Policy>1079T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos& l)1080{1081 BOOST_MATH_STD_USING1082 1083 if((a < 1) && (z < 1))1084 {1085 // No overflow possible since the power terms tend to unity as a,z -> 01086 return pow(z, a) * exp(-z) / boost::math::tgamma(a, pol);1087 }1088 else if(a > minimum_argument_for_bernoulli_recursion<T>())1089 {1090 T scaled_gamma = scaled_tgamma_no_lanczos(a, pol);1091 T power_term = pow(z / a, a / 2);1092 T a_minus_z = a - z;1093 if ((0 == power_term) || (fabs(a_minus_z) > tools::log_max_value<T>()))1094 {1095 // The result is probably zero, but we need to be sure:1096 return exp(a * log(z / a) + a_minus_z - log(scaled_gamma));1097 }1098 return (power_term * exp(a_minus_z)) * (power_term / scaled_gamma);1099 }1100 else1101 {1102 //1103 // Usual case is to calculate the prefix at a+shift and recurse down1104 // to the value we want:1105 //1106 const int min_z = minimum_argument_for_bernoulli_recursion<T>();1107 long shift = 1 + ltrunc(min_z - a);1108 T result = regularised_gamma_prefix(T(a + shift), z, pol, l);1109 if (result != 0)1110 {1111 for (long i = 0; i < shift; ++i)1112 {1113 result /= z;1114 result *= a + i;1115 }1116 return result;1117 }1118 else1119 {1120 //1121 // We failed, most probably we have z << 1, try again, this time1122 // we calculate z^a e^-z / tgamma(a+shift), combining power terms1123 // as we go. And again recurse down to the result.1124 //1125 T scaled_gamma = scaled_tgamma_no_lanczos(T(a + shift), pol);1126 T power_term_1 = pow(T(z / (a + shift)), a);1127 T power_term_2 = pow(T(a + shift), T(-shift));1128 T power_term_3 = exp(a + shift - z);1129 if ((0 == power_term_1) || (0 == power_term_2) || (0 == power_term_3) || (fabs(a + shift - z) > tools::log_max_value<T>()))1130 {1131 // We have no test case that gets here, most likely the type T1132 // has a high precision but low exponent range:1133 return exp(a * log(z) - z - boost::math::lgamma(a, pol));1134 }1135 result = power_term_1 * power_term_2 * power_term_3 / scaled_gamma;1136 for (long i = 0; i < shift; ++i)1137 {1138 result *= a + i;1139 }1140 return result;1141 }1142 }1143}1144 1145#endif // BOOST_MATH_HAS_GPU_SUPPORT1146 1147//1148// Upper gamma fraction for very small a:1149//1150template <class T, class Policy>1151BOOST_MATH_GPU_ENABLED inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)1152{1153 BOOST_MATH_STD_USING // ADL of std functions.1154 //1155 // Compute the full upper fraction (Q) when a is very small:1156 //1157 1158 #ifdef BOOST_MATH_HAS_NVRTC1159 typedef typename tools::promote_args<T>::type result_type;1160 typedef typename policies::evaluation<result_type, Policy>::type value_type;1161 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;1162 T result {detail::tgammap1m1_imp(static_cast<value_type>(a), pol, evaluation_type())};1163 #else1164 T result { boost::math::tgamma1pm1(a, pol) };1165 #endif1166 1167 if(pgam)1168 *pgam = (result + 1) / a;1169 T p = boost::math::powm1(x, a, pol);1170 result -= p;1171 result /= a;1172 detail::small_gamma2_series<T> s(a, x);1173 boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;1174 p += 1;1175 if(pderivative)1176 *pderivative = p / (*pgam * exp(x));1177 T init_value = invert ? *pgam : 0;1178 result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);1179 policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);1180 if(invert)1181 result = -result;1182 return result;1183}1184//1185// Upper gamma fraction for integer a:1186//1187template <class T, class Policy>1188BOOST_MATH_GPU_ENABLED inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)1189{1190 //1191 // Calculates normalised Q when a is an integer:1192 //1193 BOOST_MATH_STD_USING1194 T e = exp(-x);1195 T sum = e;1196 if(sum != 0)1197 {1198 T term = sum;1199 for(unsigned n = 1; n < a; ++n)1200 {1201 term /= n;1202 term *= x;1203 sum += term;1204 }1205 }1206 if(pderivative)1207 {1208 *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));1209 }1210 return sum;1211}1212//1213// Upper gamma fraction for half integer a:1214//1215template <class T, class Policy>1216BOOST_MATH_GPU_ENABLED T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)1217{1218 //1219 // Calculates normalised Q when a is a half-integer:1220 //1221 BOOST_MATH_STD_USING1222 1223 #ifdef BOOST_MATH_HAS_NVRTC1224 T e;1225 if (boost::math::is_same_v<T, float>)1226 {1227 e = ::erfcf(::sqrtf(x));1228 }1229 else1230 {1231 e = ::erfc(::sqrt(x));1232 }1233 #else1234 T e = boost::math::erfc(sqrt(x), pol);1235 #endif1236 1237 if((e != 0) && (a > 1))1238 {1239 T term = exp(-x) / sqrt(constants::pi<T>() * x);1240 term *= x;1241 static const T half = T(1) / 2; // LCOV_EXCL_LINE1242 term /= half;1243 T sum = term;1244 for(unsigned n = 2; n < a; ++n)1245 {1246 term /= n - half;1247 term *= x;1248 sum += term;1249 }1250 e += sum;1251 if(p_derivative)1252 {1253 *p_derivative = 0;1254 }1255 }1256 else if(p_derivative)1257 {1258 // We'll be dividing by x later, so calculate derivative * x:1259 *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();1260 }1261 return e;1262}1263//1264// Asymptotic approximation for large argument, see: https://dlmf.nist.gov/8.11#E21265//1266template <class T>1267struct incomplete_tgamma_large_x_series1268{1269 typedef T result_type;1270 BOOST_MATH_GPU_ENABLED incomplete_tgamma_large_x_series(const T& a, const T& x)1271 : a_poch(a - 1), z(x), term(1) {}1272 BOOST_MATH_GPU_ENABLED T operator()()1273 {1274 T result = term;1275 term *= a_poch / z;1276 a_poch -= 1;1277 return result;1278 }1279 T a_poch, z, term;1280};1281 1282template <class T, class Policy>1283BOOST_MATH_GPU_ENABLED T incomplete_tgamma_large_x(const T& a, const T& x, const Policy& pol)1284{1285 BOOST_MATH_STD_USING1286 incomplete_tgamma_large_x_series<T> s(a, x);1287 boost::math::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();1288 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);1289 boost::math::policies::check_series_iterations<T>("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol);1290 return result;1291}1292 1293 1294//1295// Main incomplete gamma entry point, handles all four incomplete gamma's:1296//1297template <class T, class Policy>1298BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, bool invert,1299 const Policy& pol, T* p_derivative)1300{1301 BOOST_MATH_STD_USING1302 1303 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;1304 1305 T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used1306 1307 BOOST_MATH_ASSERT((p_derivative == nullptr) || normalised);1308 1309 bool is_int, is_half_int;1310 bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());1311 if(is_small_a)1312 {1313 T fa = floor(a);1314 is_int = (fa == a);1315 is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);1316 }1317 else1318 {1319 is_int = is_half_int = false;1320 }1321 1322 int eval_method;1323 1324 if (x == 0)1325 {1326 eval_method = 2;1327 }1328 else if(is_int && (x > 0.6))1329 {1330 // calculate Q via finite sum:1331 invert = !invert;1332 eval_method = 0;1333 }1334 else if(is_half_int && (x > 0.2))1335 {1336 // calculate Q via finite sum for half integer a:1337 invert = !invert;1338 eval_method = 1;1339 }1340 else if((x < tools::root_epsilon<T>()) && (a > 1))1341 {1342 eval_method = 6;1343 }1344 else if ((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1)))1345 {1346 // calculate Q via asymptotic approximation:1347 invert = !invert;1348 eval_method = 7;1349 }1350 else if(x < T(0.5))1351 {1352 //1353 // Changeover criterion chosen to give a changeover at Q ~ 0.331354 //1355 if(T(-0.4) / log(x) < a)1356 {1357 eval_method = 2;1358 }1359 else1360 {1361 eval_method = 3;1362 }1363 }1364 else if(x < T(1.1))1365 {1366 //1367 // Changeover here occurs when P ~ 0.75 or Q ~ 0.25:1368 //1369 if(x * 0.75f < a)1370 {1371 eval_method = 2;1372 }1373 else1374 {1375 eval_method = 3;1376 }1377 }1378 else1379 {1380 //1381 // Begin by testing whether we're in the "bad" zone1382 // where the result will be near 0.5 and the usual1383 // series and continued fractions are slow to converge:1384 //1385 bool use_temme = false;1386 if(normalised && boost::math::numeric_limits<T>::is_specialized && (a > 20))1387 {1388 T sigma = fabs((x-a)/a);1389 if((a > 200) && (policies::digits<T, Policy>() <= 113))1390 {1391 //1392 // This limit is chosen so that we use Temme's expansion1393 // only if the result would be larger than about 10^-6.1394 // Below that the regular series and continued fractions1395 // converge OK, and if we use Temme's method we get increasing1396 // errors from the dominant erfc term as it's (inexact) argument1397 // increases in magnitude.1398 //1399 if(20 / a > sigma * sigma)1400 use_temme = true;1401 }1402 else if(policies::digits<T, Policy>() <= 64)1403 {1404 // Note in this zone we can't use Temme's expansion for1405 // types longer than an 80-bit real:1406 // it would require too many terms in the polynomials.1407 if(sigma < 0.4)1408 use_temme = true;1409 }1410 }1411 if(use_temme)1412 {1413 eval_method = 5;1414 }1415 else1416 {1417 //1418 // Regular case where the result will not be too close to 0.5.1419 //1420 // Changeover here occurs at P ~ Q ~ 0.51421 // Note that series computation of P is about x2 faster than continued fraction1422 // calculation of Q, so try and use the CF only when really necessary, especially1423 // for small x.1424 //1425 if(x - (1 / (3 * x)) < a)1426 {1427 eval_method = 2;1428 }1429 else1430 {1431 eval_method = 4;1432 invert = !invert;1433 }1434 }1435 }1436 1437 switch(eval_method)1438 {1439 case 0:1440 {1441 result = finite_gamma_q(a, x, pol, p_derivative);1442 if(!normalised)1443 {1444 #ifdef BOOST_MATH_HAS_NVRTC1445 if (boost::math::is_same_v<T, float>)1446 {1447 result *= ::tgammaf(a);1448 }1449 else1450 {1451 result *= ::tgamma(a);1452 }1453 #else1454 result *= boost::math::tgamma(a, pol);1455 #endif1456 }1457 break;1458 }1459 case 1:1460 {1461 result = finite_half_gamma_q(a, x, p_derivative, pol);1462 if(!normalised)1463 {1464 #ifdef BOOST_MATH_HAS_NVRTC1465 if (boost::math::is_same_v<T, float>)1466 {1467 result *= ::tgammaf(a);1468 }1469 else1470 {1471 result *= ::tgamma(a);1472 }1473 #else1474 result *= boost::math::tgamma(a, pol);1475 #endif1476 }1477 if(p_derivative && (*p_derivative == 0))1478 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());1479 break;1480 }1481 case 2:1482 {1483 // Compute P:1484 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);1485 if(p_derivative)1486 *p_derivative = result;1487 if(result != 0)1488 {1489 //1490 // If we're going to be inverting the result then we can1491 // reduce the number of series evaluations by quite1492 // a few iterations if we set an initial value for the1493 // series sum based on what we'll end up subtracting it from1494 // at the end.1495 // Have to be careful though that this optimization doesn't1496 // lead to spurious numeric overflow. Note that the1497 // scary/expensive overflow checks below are more often1498 // than not bypassed in practice for "sensible" input1499 // values:1500 //1501 T init_value = 0;1502 bool optimised_invert = false;1503 if(invert)1504 {1505 #ifdef BOOST_MATH_HAS_NVRTC1506 if (boost::math::is_same_v<T, float>)1507 {1508 init_value = (normalised ? T(1) : ::tgammaf(a));1509 }1510 else1511 {1512 init_value = (normalised ? T(1) : ::tgamma(a));1513 }1514 #else1515 init_value = (normalised ? T(1) : boost::math::tgamma(a, pol));1516 #endif1517 1518 if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))1519 {1520 init_value /= result;1521 if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value))1522 {1523 init_value *= -a;1524 optimised_invert = true;1525 }1526 else1527 init_value = 0; // LCOV_EXCL_LINE Unreachable for any "sensible" floating point type.1528 }1529 else1530 init_value = 0;1531 }1532 result *= detail::lower_gamma_series(a, x, pol, init_value) / a;1533 if(optimised_invert)1534 {1535 invert = false;1536 result = -result;1537 }1538 }1539 break;1540 }1541 case 3:1542 {1543 // Compute Q:1544 invert = !invert;1545 T g{};1546 result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);1547 invert = false;1548 if(normalised)1549 result /= g;1550 break;1551 }1552 case 4:1553 {1554 // Compute Q:1555 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);1556 if(p_derivative)1557 *p_derivative = result;1558 if(result != 0)1559 result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());1560 break;1561 }1562 case 5:1563 {1564 //1565 // Use compile time dispatch to the appropriate1566 // Temme asymptotic expansion. This may be dead code1567 // if T does not have numeric limits support, or has1568 // too many digits for the most precise version of1569 // these expansions, in that case we'll be calling1570 // an empty function.1571 //1572 typedef typename policies::precision<T, Policy>::type precision_type;1573 1574 typedef boost::math::integral_constant<int,1575 precision_type::value <= 0 ? 0 :1576 precision_type::value <= 53 ? 53 :1577 precision_type::value <= 64 ? 64 :1578 precision_type::value <= 113 ? 113 : 01579 > tag_type;1580 1581 result = igamma_temme_large(a, x, pol, tag_type());1582 if(x >= a)1583 invert = !invert;1584 if(p_derivative)1585 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());1586 break;1587 }1588 case 6:1589 {1590 // x is so small that P is necessarily very small too,1591 // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/1592 if(!normalised)1593 result = pow(x, a) / (a);1594 else1595 {1596#ifndef BOOST_MATH_NO_EXCEPTIONS1597 try1598 {1599#endif1600 #ifdef BOOST_MATH_HAS_NVRTC1601 if (boost::math::is_same_v<T, float>)1602 {1603 result = ::powf(x, a) / ::tgammaf(a + 1);1604 }1605 else1606 {1607 result = ::pow(x, a) / ::tgamma(a + 1);1608 }1609 #else1610 result = pow(x, a) / boost::math::tgamma(a + 1, pol);1611 #endif1612#ifndef BOOST_MATH_NO_EXCEPTIONS1613 }1614 catch (const std::overflow_error&)1615 {1616 result = 0;1617 }1618#endif1619 }1620 result *= 1 - a * x / (a + 1);1621 if (p_derivative)1622 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());1623 break;1624 }1625 case 7:1626 {1627 // x is large,1628 // Compute Q:1629 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);1630 if (p_derivative)1631 *p_derivative = result;1632 result /= x;1633 if (result != 0)1634 result *= incomplete_tgamma_large_x(a, x, pol);1635 break;1636 }1637 }1638 1639 if(normalised && (result > 1))1640 result = 1;1641 if(invert)1642 {1643 #ifdef BOOST_MATH_HAS_NVRTC1644 T gam;1645 if (boost::math::is_same_v<T, float>)1646 {1647 gam = normalised ? T(1) : ::tgammaf(a);1648 }1649 else1650 {1651 gam = normalised ? T(1) : ::tgamma(a);1652 }1653 #else1654 T gam = normalised ? T(1) : boost::math::tgamma(a, pol);1655 #endif1656 result = gam - result;1657 }1658 if(p_derivative)1659 {1660 if((x == 0) || ((x < 1) && (tools::max_value<T>() * x < *p_derivative)))1661 {1662 // overflow, just return an arbitrarily large value:1663 *p_derivative = tools::max_value<T>() / 2;1664 }1665 else1666 *p_derivative /= x;1667 }1668 1669 return result;1670}1671 1672// Need to implement this dispatch to avoid recursion for device compilers1673template <class T, class Policy>1674BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,1675 const Policy& pol, T* p_derivative)1676{1677 constexpr auto function = "boost::math::gamma_p<%1%>(%1%, %1%)";1678 if(a <= 0)1679 return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);1680 if(x < 0)1681 return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);1682 1683 BOOST_MATH_STD_USING1684 1685 1686 T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used1687 1688 if(x > 0 && a >= max_factorial<T>::value && !normalised)1689 {1690 //1691 // When we're computing the non-normalized incomplete gamma1692 // and a is large the result is rather hard to compute unless1693 // we use logs. There are really two options - if x is a long1694 // way from a in value then we can reliably use methods 2 and 41695 // below in logarithmic form and go straight to the result.1696 // Otherwise we let the regularized gamma take the strain1697 // (the result is unlikely to underflow in the central region anyway)1698 // and combine with lgamma in the hopes that we get a finite result.1699 //1700 if(invert && (a * 4 < x))1701 {1702 // This is method 4 below, done in logs:1703 result = a * log(x) - x;1704 BOOST_MATH_ASSERT(p_derivative == nullptr);1705 // Not currently used for non-normalized igamma:1706 //if(p_derivative)1707 // *p_derivative = exp(result);1708 result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));1709 }1710 else if(!invert && (a > 4 * x))1711 {1712 // This is method 2 below, done in logs:1713 result = a * log(x) - x;1714 BOOST_MATH_ASSERT(p_derivative == nullptr);1715 // Not currently used for non-normalized igamma:1716 //if(p_derivative)1717 // *p_derivative = exp(result);1718 T init_value = 0;1719 result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);1720 }1721 else1722 {1723 result = gamma_incomplete_imp_final(T(a), T(x), true, invert, pol, p_derivative);1724 if(result == 0)1725 {1726 if(invert)1727 {1728 // Try http://functions.wolfram.com/06.06.06.0039.011729 result = 1 + 1 / (12 * a) + 1 / (288 * a * a);1730 result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>());1731 BOOST_MATH_ASSERT(p_derivative == nullptr);1732 // Not currently used for non-normalized igamma:1733 //if(p_derivative)1734 // *p_derivative = exp(a * log(x) - x);1735 }1736 else1737 {1738 // This is method 2 below, done in logs, we're really outside the1739 // range of this method, but since the result is almost certainly1740 // infinite, we should probably be OK:1741 result = a * log(x) - x;1742 BOOST_MATH_ASSERT(p_derivative == nullptr);1743 // Not currently used for non-normalized igamma:1744 //if(p_derivative)1745 // *p_derivative = exp(result);1746 T init_value = 0;1747 result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);1748 }1749 }1750 else1751 {1752 #ifdef BOOST_MATH_HAS_NVRTC1753 if (boost::math::is_same_v<T, float>)1754 {1755 result = ::logf(result) + ::lgammaf(a);1756 }1757 else1758 {1759 result = ::log(result) + ::lgamma(a);1760 }1761 #else1762 result = log(result) + boost::math::lgamma(a, pol);1763 #endif1764 }1765 }1766 if(result > tools::log_max_value<T>())1767 return policies::raise_overflow_error<T>(function, nullptr, pol);1768 return exp(result);1769 }1770 1771 // If no special handling is required then we proceeds as normal1772 return gamma_incomplete_imp_final(T(a), T(x), normalised, invert, pol, p_derivative);1773}1774 1775//1776// Ratios of two gamma functions:1777//1778template <class T, class Policy, class Lanczos>1779BOOST_MATH_GPU_ENABLED T tgamma_delta_ratio_imp_lanczos_final(T z, T delta, const Policy& pol, const Lanczos&)1780{1781 BOOST_MATH_STD_USING1782 1783 T zgh = static_cast<T>(z + T(Lanczos::g()) - constants::half<T>());1784 T result{};1785 if(z + delta == z)1786 {1787 // Given delta < z * eps1788 // and zgh > z1789 // Then this must follow:1790 BOOST_MATH_ASSERT(fabs(delta / zgh) < boost::math::tools::epsilon<T>());1791 // We have:1792 // result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));1793 // 0.5 - z == -z1794 // log1p(delta / zgh) = delta / zgh = delta / z1795 // multiplying we get -delta.1796 result = exp(-delta);1797 }1798 else1799 {1800 if(fabs(delta) < 10)1801 {1802 result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));1803 }1804 else1805 {1806 result = pow(T(zgh / (zgh + delta)), T(z - constants::half<T>()));1807 }1808 // Split the calculation up to avoid spurious overflow:1809 result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));1810 }1811 result *= pow(T(constants::e<T>() / (zgh + delta)), delta);1812 return result;1813}1814 1815template <class T, class Policy, class Lanczos>1816BOOST_MATH_GPU_ENABLED T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l)1817{1818 BOOST_MATH_STD_USING1819 1820 if(z < tools::epsilon<T>())1821 {1822 //1823 // We get spurious numeric overflow unless we're very careful, this1824 // can occur either inside Lanczos::lanczos_sum(z) or in the1825 // final combination of terms, to avoid this, split the product up1826 // into 2 (or 3) parts:1827 //1828 // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta1829 // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial1830 //1831 if(boost::math::max_factorial<T>::value < delta)1832 {1833 T ratio = tgamma_delta_ratio_imp_lanczos_final(T(delta), T(boost::math::max_factorial<T>::value - delta), pol, l);1834 ratio *= z;1835 ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1);1836 return 1 / ratio;1837 }1838 else1839 {1840 #ifdef BOOST_MATH_HAS_NVRTC1841 if (boost::math::is_same_v<T, float>)1842 {1843 return 1 / (z * ::tgammaf(z + delta));1844 }1845 else1846 {1847 return 1 / (z * ::tgamma(z + delta));1848 }1849 #else1850 return 1 / (z * boost::math::tgamma(z + delta, pol));1851 #endif1852 }1853 }1854 1855 return tgamma_delta_ratio_imp_lanczos_final(T(z), T(delta), pol, l);1856}1857 1858//1859// And again without Lanczos support this time:1860//1861#ifndef BOOST_MATH_HAS_GPU_SUPPORT1862 1863template <class T, class Policy>1864T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos& l)1865{1866 BOOST_MATH_STD_USING1867 1868 //1869 // We adjust z and delta so that both z and z+delta are large enough for1870 // Sterling's approximation to hold. We can then calculate the ratio1871 // for the adjusted values, and rescale back down to z and z+delta.1872 //1873 // Get the required shifts first:1874 //1875 long numerator_shift = 0;1876 long denominator_shift = 0;1877 const int min_z = minimum_argument_for_bernoulli_recursion<T>();1878 1879 if (min_z > z)1880 numerator_shift = 1 + ltrunc(min_z - z);1881 if (min_z > z + delta)1882 denominator_shift = 1 + ltrunc(min_z - z - delta);1883 //1884 // If the shifts are zero, then we can just combine scaled tgamma's1885 // and combine the remaining terms:1886 //1887 if (numerator_shift == 0 && denominator_shift == 0)1888 {1889 T scaled_tgamma_num = scaled_tgamma_no_lanczos(z, pol);1890 T scaled_tgamma_denom = scaled_tgamma_no_lanczos(T(z + delta), pol);1891 T result = scaled_tgamma_num / scaled_tgamma_denom;1892 result *= exp(z * boost::math::log1p(-delta / (z + delta), pol)) * pow(T((delta + z) / constants::e<T>()), -delta);1893 return result;1894 }1895 //1896 // We're going to have to rescale first, get the adjusted z and delta values,1897 // plus the ratio for the adjusted values:1898 //1899 T zz = z + numerator_shift;1900 T dd = delta - (numerator_shift - denominator_shift);1901 T ratio = tgamma_delta_ratio_imp_lanczos(zz, dd, pol, l);1902 //1903 // Use gamma recurrence relations to get back to the original1904 // z and z+delta:1905 //1906 for (long long i = 0; i < numerator_shift; ++i)1907 {1908 ratio /= (z + i);1909 if (i < denominator_shift)1910 ratio *= (z + delta + i);1911 }1912 for (long long i = numerator_shift; i < denominator_shift; ++i)1913 {1914 ratio *= (z + delta + i);1915 }1916 return ratio;1917}1918 1919#endif1920 1921template <class T, class Policy>1922BOOST_MATH_GPU_ENABLED T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)1923{1924 BOOST_MATH_STD_USING1925 1926 if((z <= 0) || (z + delta <= 0))1927 {1928 // This isn't very sophisticated, or accurate, but it does work:1929 #ifdef BOOST_MATH_HAS_NVRTC1930 if (boost::math::is_same_v<T, float>)1931 {1932 return ::tgammaf(z) / ::tgammaf(z + delta);1933 }1934 else1935 {1936 return ::tgamma(z) / ::tgamma(z + delta);1937 }1938 #else1939 return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol);1940 #endif1941 }1942 1943 if(floor(delta) == delta)1944 {1945 if(floor(z) == z)1946 {1947 //1948 // Both z and delta are integers, see if we can just use table lookup1949 // of the factorials to get the result:1950 //1951 if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))1952 {1953 return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);1954 }1955 }1956 if(fabs(delta) < 20)1957 {1958 //1959 // delta is a small integer, we can use a finite product:1960 //1961 if(delta == 0)1962 return 1;1963 if(delta < 0)1964 {1965 z -= 1;1966 T result = z;1967 while(0 != (delta += 1))1968 {1969 z -= 1;1970 result *= z;1971 }1972 return result;1973 }1974 else1975 {1976 T result = 1 / z;1977 while(0 != (delta -= 1))1978 {1979 z += 1;1980 result /= z;1981 }1982 return result;1983 }1984 }1985 }1986 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;1987 return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());1988}1989 1990template <class T, class Policy>1991BOOST_MATH_GPU_ENABLED T tgamma_ratio_imp(T x, T y, const Policy& pol)1992{1993 BOOST_MATH_STD_USING1994 1995 if((x <= 0) || (boost::math::isinf)(x))1996 return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", x, pol);1997 if((y <= 0) || (boost::math::isinf)(y))1998 return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol);1999 2000 // We don't need to worry about the denorm case on device2001 // And this has the added bonus of removing recursion2002 #ifndef BOOST_MATH_HAS_GPU_SUPPORT2003 if(x <= tools::min_value<T>())2004 {2005 // Special case for denorms...Ugh.2006 T shift = ldexp(T(1), tools::digits<T>());2007 return shift * tgamma_ratio_imp(T(x * shift), y, pol);2008 }2009 #endif2010 2011 if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))2012 {2013 // Rather than subtracting values, lets just call the gamma functions directly:2014 #ifdef BOOST_MATH_HAS_NVRTC2015 if (boost::math::is_same_v<T, float>)2016 {2017 return ::tgammaf(x) / ::tgammaf(y);2018 }2019 else2020 {2021 return ::tgamma(x) / ::tgamma(y);2022 }2023 #else2024 return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);2025 #endif2026 }2027 T prefix = 1;2028 if(x < 1)2029 {2030 if(y < 2 * max_factorial<T>::value)2031 {2032 // We need to sidestep on x as well, otherwise we'll underflow2033 // before we get to factor in the prefix term:2034 prefix /= x;2035 x += 1;2036 while(y >= max_factorial<T>::value)2037 {2038 y -= 1;2039 prefix /= y;2040 }2041 2042 #ifdef BOOST_MATH_HAS_NVRTC2043 if (boost::math::is_same_v<T, float>)2044 {2045 return prefix * ::tgammaf(x) / ::tgammaf(y);2046 }2047 else2048 {2049 return prefix * ::tgamma(x) / ::tgamma(y);2050 }2051 #else2052 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);2053 #endif2054 }2055 //2056 // result is almost certainly going to underflow to zero, try logs just in case:2057 //2058 #ifdef BOOST_MATH_HAS_NVRTC2059 if (boost::math::is_same_v<T, float>)2060 {2061 return ::expf(::lgammaf(x) - ::lgammaf(y));2062 }2063 else2064 {2065 return ::exp(::lgamma(x) - ::lgamma(y));2066 }2067 #else2068 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));2069 #endif2070 }2071 if(y < 1)2072 {2073 if(x < 2 * max_factorial<T>::value)2074 {2075 // We need to sidestep on y as well, otherwise we'll overflow2076 // before we get to factor in the prefix term:2077 prefix *= y;2078 y += 1;2079 while(x >= max_factorial<T>::value)2080 {2081 x -= 1;2082 prefix *= x;2083 }2084 2085 #ifdef BOOST_MATH_HAS_NVRTC2086 if (boost::math::is_same_v<T, float>)2087 {2088 return prefix * ::tgammaf(x) / ::tgammaf(y);2089 }2090 else2091 {2092 return prefix * ::tgamma(x) / ::tgamma(y);2093 }2094 #else2095 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);2096 #endif2097 }2098 //2099 // Result will almost certainly overflow, try logs just in case:2100 //2101 BOOST_MATH_IF_CONSTEXPR(boost::math::is_same<T, float>::value || boost::math::is_same<T, double>::value)2102 {2103 // straight to the scene of the accident, since the result is larger than max_factorial:2104 return policies::raise_overflow_error<T>("tgamma_ratio", nullptr, pol);2105 }2106 else2107 {2108 #ifdef BOOST_MATH_HAS_NVRTC2109 if (boost::math::is_same_v<T, float>)2110 {2111 prefix = ::lgammaf(x) - ::lgammaf(y);2112 }2113 else2114 {2115 prefix = ::lgamma(x) - ::lgamma(y);2116 }2117 #else2118 prefix = boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol);2119 #endif2120 if (prefix > boost::math::tools::log_max_value<T>())2121 return policies::raise_overflow_error<T>("tgamma_ratio", nullptr, pol);2122 //2123 // This is unreachable, unless max_factorial is small compared to the exponent2124 // range of the type, ie multiprecision types only here...2125 //2126 return exp(prefix); // LCOV_EXCL_LINE2127 }2128 }2129 //2130 // Regular case, x and y both large and similar in magnitude:2131 //2132 #ifdef BOOST_MATH_HAS_NVRTC2133 return detail::tgamma_delta_ratio_imp(x, y - x, pol);2134 #else2135 return boost::math::tgamma_delta_ratio(x, y - x, pol);2136 #endif2137}2138 2139template <class T, class Policy>2140BOOST_MATH_GPU_ENABLED T gamma_p_derivative_imp(T a, T x, const Policy& pol)2141{2142 BOOST_MATH_STD_USING2143 //2144 // Usual error checks first:2145 //2146 if(a <= 0)2147 return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);2148 if(x < 0)2149 return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);2150 //2151 // Now special cases:2152 //2153 if(x == 0)2154 {2155 return (a > 1) ? T(0) :2156 (a == 1) ? T(1) : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", nullptr, pol);2157 }2158 //2159 // Normal case:2160 //2161 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;2162 T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());2163 /*2164 * Derivative goes to zero as x -> 0, this should be unreachable:2165 * 2166 if((x < 1) && (tools::max_value<T>() * x < f1))2167 {2168 // overflow:2169 return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", nullptr, pol);2170 }2171 */2172 if(f1 == 0)2173 {2174 // Underflow in calculation, use logs instead:2175 #ifdef BOOST_MATH_HAS_NVRTC2176 if (boost::math::is_same_v<T, float>)2177 {2178 f1 = a * ::logf(x) - x - ::lgammaf(a) - ::logf(x);2179 }2180 else2181 {2182 f1 = a * ::log(x) - x - ::lgamma(a) - ::log(x);2183 }2184 #else2185 f1 = a * log(x) - x - lgamma(a, pol) - log(x);2186 #endif2187 f1 = exp(f1);2188 }2189 else2190 f1 /= x;2191 2192 return f1;2193}2194 2195template <class T, class Policy>2196BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type2197 tgamma(T z, const Policy& /* pol */, const boost::math::true_type)2198{2199 BOOST_FPU_EXCEPTION_GUARD2200 typedef typename tools::promote_args<T>::type result_type;2201 typedef typename policies::evaluation<result_type, Policy>::type value_type;2202 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;2203 typedef typename policies::normalise<2204 Policy,2205 policies::promote_float<false>,2206 policies::promote_double<false>,2207 policies::discrete_quantile<>,2208 policies::assert_undefined<> >::type forwarding_policy;2209 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)");2210}2211 2212template <class T1, class T2, class Policy>2213BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2214 tgamma(T1 a, T2 z, const Policy&, const boost::math::false_type)2215{2216 BOOST_FPU_EXCEPTION_GUARD2217 typedef tools::promote_args_t<T1, T2> result_type;2218 typedef typename policies::evaluation<result_type, Policy>::type value_type;2219 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;2220 typedef typename policies::normalise<2221 Policy,2222 policies::promote_float<false>,2223 policies::promote_double<false>,2224 policies::discrete_quantile<>,2225 policies::assert_undefined<> >::type forwarding_policy;2226 2227 return policies::checked_narrowing_cast<result_type, forwarding_policy>(2228 detail::gamma_incomplete_imp(static_cast<value_type>(a),2229 static_cast<value_type>(z), false, true,2230 forwarding_policy(), static_cast<value_type*>(nullptr)), "boost::math::tgamma<%1%>(%1%, %1%)");2231}2232 2233template <class T1, class T2>2234BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2235 tgamma(T1 a, T2 z, const boost::math::false_type& tag)2236{2237 return tgamma(a, z, policies::policy<>(), tag);2238}2239 2240 2241} // namespace detail2242 2243template <class T, class Policy>2244BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type2245 lgamma(T z, int* sign, const Policy&)2246{2247 BOOST_FPU_EXCEPTION_GUARD2248 typedef typename tools::promote_args<T>::type result_type;2249 typedef typename policies::evaluation<result_type, Policy>::type value_type;2250 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;2251 typedef typename policies::normalise<2252 Policy,2253 policies::promote_float<false>,2254 policies::promote_double<false>,2255 policies::discrete_quantile<>,2256 policies::assert_undefined<> >::type forwarding_policy;2257 2258 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lgamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)");2259}2260 2261template <class T>2262BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type2263 lgamma(T z, int* sign)2264{2265 return lgamma(z, sign, policies::policy<>());2266}2267 2268template <class T, class Policy>2269BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type2270 lgamma(T x, const Policy& pol)2271{2272 return ::boost::math::lgamma(x, nullptr, pol);2273}2274 2275template <class T>2276BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type2277 lgamma(T x)2278{2279 return ::boost::math::lgamma(x, nullptr, policies::policy<>());2280}2281 2282template <class T, class Policy>2283BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type2284 tgamma1pm1(T z, const Policy& /* pol */)2285{2286 BOOST_FPU_EXCEPTION_GUARD2287 typedef typename tools::promote_args<T>::type result_type;2288 typedef typename policies::evaluation<result_type, Policy>::type value_type;2289 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;2290 typedef typename policies::normalise<2291 Policy,2292 policies::promote_float<false>,2293 policies::promote_double<false>,2294 policies::discrete_quantile<>,2295 policies::assert_undefined<> >::type forwarding_policy;2296 2297 return policies::checked_narrowing_cast<typename boost::math::remove_cv<result_type>::type, forwarding_policy>(detail::tgammap1m1_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)");2298}2299 2300template <class T>2301BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type2302 tgamma1pm1(T z)2303{2304 return tgamma1pm1(z, policies::policy<>());2305}2306 2307//2308// Full upper incomplete gamma:2309//2310template <class T1, class T2>2311BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2312 tgamma(T1 a, T2 z)2313{2314 //2315 // Type T2 could be a policy object, or a value, select the2316 // right overload based on T2:2317 //2318 using maybe_policy = typename policies::is_policy<T2>::type;2319 using result_type = tools::promote_args_t<T1, T2>;2320 return static_cast<result_type>(detail::tgamma(a, z, maybe_policy()));2321}2322template <class T1, class T2, class Policy>2323BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2324 tgamma(T1 a, T2 z, const Policy& pol)2325{2326 using result_type = tools::promote_args_t<T1, T2>;2327 return static_cast<result_type>(detail::tgamma(a, z, pol, boost::math::false_type()));2328}2329template <class T>2330BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type2331 tgamma(T z)2332{2333 return tgamma(z, policies::policy<>());2334}2335//2336// Full lower incomplete gamma:2337//2338template <class T1, class T2, class Policy>2339BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2340 tgamma_lower(T1 a, T2 z, const Policy&)2341{2342 BOOST_FPU_EXCEPTION_GUARD2343 typedef tools::promote_args_t<T1, T2> result_type;2344 typedef typename policies::evaluation<result_type, Policy>::type value_type;2345 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;2346 typedef typename policies::normalise<2347 Policy,2348 policies::promote_float<false>,2349 policies::promote_double<false>,2350 policies::discrete_quantile<>,2351 policies::assert_undefined<> >::type forwarding_policy;2352 2353 return policies::checked_narrowing_cast<result_type, forwarding_policy>(2354 detail::gamma_incomplete_imp(static_cast<value_type>(a),2355 static_cast<value_type>(z), false, false,2356 forwarding_policy(), static_cast<value_type*>(nullptr)), "tgamma_lower<%1%>(%1%, %1%)");2357}2358template <class T1, class T2>2359BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2360 tgamma_lower(T1 a, T2 z)2361{2362 return tgamma_lower(a, z, policies::policy<>());2363}2364//2365// Regularised upper incomplete gamma:2366//2367template <class T1, class T2, class Policy>2368BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2369 gamma_q(T1 a, T2 z, const Policy& /* pol */)2370{2371 BOOST_FPU_EXCEPTION_GUARD2372 typedef tools::promote_args_t<T1, T2> result_type;2373 typedef typename policies::evaluation<result_type, Policy>::type value_type;2374 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;2375 typedef typename policies::normalise<2376 Policy,2377 policies::promote_float<false>,2378 policies::promote_double<false>,2379 policies::discrete_quantile<>,2380 policies::assert_undefined<> >::type forwarding_policy;2381 2382 return policies::checked_narrowing_cast<result_type, forwarding_policy>(2383 detail::gamma_incomplete_imp(static_cast<value_type>(a),2384 static_cast<value_type>(z), true, true,2385 forwarding_policy(), static_cast<value_type*>(nullptr)), "gamma_q<%1%>(%1%, %1%)");2386}2387template <class T1, class T2>2388BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2389 gamma_q(T1 a, T2 z)2390{2391 return gamma_q(a, z, policies::policy<>());2392}2393//2394// Regularised lower incomplete gamma:2395//2396template <class T1, class T2, class Policy>2397BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2398 gamma_p(T1 a, T2 z, const Policy&)2399{2400 BOOST_FPU_EXCEPTION_GUARD2401 typedef tools::promote_args_t<T1, T2> result_type;2402 typedef typename policies::evaluation<result_type, Policy>::type value_type;2403 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;2404 typedef typename policies::normalise<2405 Policy,2406 policies::promote_float<false>,2407 policies::promote_double<false>,2408 policies::discrete_quantile<>,2409 policies::assert_undefined<> >::type forwarding_policy;2410 2411 return policies::checked_narrowing_cast<result_type, forwarding_policy>(2412 detail::gamma_incomplete_imp(static_cast<value_type>(a),2413 static_cast<value_type>(z), true, false,2414 forwarding_policy(), static_cast<value_type*>(nullptr)), "gamma_p<%1%>(%1%, %1%)");2415}2416template <class T1, class T2>2417BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2418 gamma_p(T1 a, T2 z)2419{2420 return gamma_p(a, z, policies::policy<>());2421}2422 2423// ratios of gamma functions:2424template <class T1, class T2, class Policy>2425BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2426 tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)2427{2428 BOOST_FPU_EXCEPTION_GUARD2429 typedef tools::promote_args_t<T1, T2> result_type;2430 typedef typename policies::evaluation<result_type, Policy>::type value_type;2431 typedef typename policies::normalise<2432 Policy,2433 policies::promote_float<false>,2434 policies::promote_double<false>,2435 policies::discrete_quantile<>,2436 policies::assert_undefined<> >::type forwarding_policy;2437 2438 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(z), static_cast<value_type>(delta), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");2439}2440template <class T1, class T2>2441BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2442 tgamma_delta_ratio(T1 z, T2 delta)2443{2444 return tgamma_delta_ratio(z, delta, policies::policy<>());2445}2446template <class T1, class T2, class Policy>2447BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2448 tgamma_ratio(T1 a, T2 b, const Policy&)2449{2450 typedef tools::promote_args_t<T1, T2> result_type;2451 typedef typename policies::evaluation<result_type, Policy>::type value_type;2452 typedef typename policies::normalise<2453 Policy,2454 policies::promote_float<false>,2455 policies::promote_double<false>,2456 policies::discrete_quantile<>,2457 policies::assert_undefined<> >::type forwarding_policy;2458 2459 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");2460}2461template <class T1, class T2>2462BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2463 tgamma_ratio(T1 a, T2 b)2464{2465 return tgamma_ratio(a, b, policies::policy<>());2466}2467 2468template <class T1, class T2, class Policy>2469BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2470 gamma_p_derivative(T1 a, T2 x, const Policy&)2471{2472 BOOST_FPU_EXCEPTION_GUARD2473 typedef tools::promote_args_t<T1, T2> result_type;2474 typedef typename policies::evaluation<result_type, Policy>::type value_type;2475 typedef typename policies::normalise<2476 Policy,2477 policies::promote_float<false>,2478 policies::promote_double<false>,2479 policies::discrete_quantile<>,2480 policies::assert_undefined<> >::type forwarding_policy;2481 2482 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_p_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(x), forwarding_policy()), "boost::math::gamma_p_derivative<%1%>(%1%, %1%)");2483}2484template <class T1, class T2>2485BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>2486 gamma_p_derivative(T1 a, T2 x)2487{2488 return gamma_p_derivative(a, x, policies::policy<>());2489}2490 2491} // namespace math2492} // namespace boost2493 2494#ifdef _MSC_VER2495# pragma warning(pop)2496#endif2497 2498#include <boost/math/special_functions/detail/igamma_inverse.hpp>2499#include <boost/math/special_functions/detail/gamma_inva.hpp>2500#include <boost/math/special_functions/erf.hpp>2501 2502#endif // BOOST_MATH_SF_GAMMA_HPP2503