1[section:igamma Incomplete Gamma Functions] 2 3[h4 Synopsis] 4 5`` 6#include <boost/math/special_functions/gamma.hpp> 7`` 8 9 namespace boost{ namespace math{ 10 11 template <class T1, class T2> 12 ``__sf_result`` gamma_p(T1 a, T2 z); 13 14 template <class T1, class T2, class ``__Policy``> 15 ``__sf_result`` gamma_p(T1 a, T2 z, const ``__Policy``&); 16 17 template <class T1, class T2> 18 ``__sf_result`` gamma_q(T1 a, T2 z); 19 20 template <class T1, class T2, class ``__Policy``> 21 ``__sf_result`` gamma_q(T1 a, T2 z, const ``__Policy``&); 22 23 template <class T1, class T2> 24 ``__sf_result`` tgamma_lower(T1 a, T2 z); 25 26 template <class T1, class T2, class ``__Policy``> 27 ``__sf_result`` tgamma_lower(T1 a, T2 z, const ``__Policy``&); 28 29 template <class T1, class T2> 30 ``__sf_result`` tgamma(T1 a, T2 z); 31 32 template <class T1, class T2, class ``__Policy``> 33 ``__sf_result`` tgamma(T1 a, T2 z, const ``__Policy``&); 34 35 }} // namespaces 36 37[h4 Description] 38 39There are four [@http://mathworld.wolfram.com/IncompleteGammaFunction.html 40incomplete gamma functions]: 41two are normalised versions (also known as /regularized/ incomplete gamma functions) 42that return values in the range [0, 1], and two are non-normalised and 43return values in the range [0, [Gamma](a)]. Users interested in statistical 44applications should use the 45[@http://mathworld.wolfram.com/RegularizedGammaFunction.html normalised versions (`gamma_p` and `gamma_q`)]. 46 47All of these functions require /a > 0/ and /z >= 0/, otherwise they return 48the result of __domain_error. 49 50[optional_policy] 51 52The return type of these functions is computed using the __arg_promotion_rules 53when T1 and T2 are different types, otherwise the return type is simply T1. 54 55 template <class T1, class T2> 56 ``__sf_result`` gamma_p(T1 a, T2 z); 57 58 template <class T1, class T2, class Policy> 59 ``__sf_result`` gamma_p(T1 a, T2 z, const ``__Policy``&); 60 61Returns the normalised lower incomplete gamma function of a and z: 62 63[equation igamma4] 64 65This function changes rapidly from 0 to 1 around the point z == a: 66 67[graph gamma_p] 68 69 template <class T1, class T2> 70 ``__sf_result`` gamma_q(T1 a, T2 z); 71 72 template <class T1, class T2, class ``__Policy``> 73 ``__sf_result`` gamma_q(T1 a, T2 z, const ``__Policy``&); 74 75Returns the normalised upper incomplete gamma function of a and z: 76 77[equation igamma3] 78 79This function changes rapidly from 1 to 0 around the point z == a: 80 81[graph gamma_q] 82 83 template <class T1, class T2> 84 ``__sf_result`` tgamma_lower(T1 a, T2 z); 85 86 template <class T1, class T2, class ``__Policy``> 87 ``__sf_result`` tgamma_lower(T1 a, T2 z, const ``__Policy``&); 88 89Returns the full (non-normalised) lower incomplete gamma function of a and z: 90 91[equation igamma2] 92 93 template <class T1, class T2> 94 ``__sf_result`` tgamma(T1 a, T2 z); 95 96 template <class T1, class T2, class ``__Policy``> 97 ``__sf_result`` tgamma(T1 a, T2 z, const ``__Policy``&); 98 99Returns the full (non-normalised) upper incomplete gamma function of a and z: 100 101[equation igamma1] 102 103[h4 Accuracy] 104 105The following tables give peak and mean relative errors in over various domains of 106a and z, along with comparisons to the __gsl and __cephes libraries. 107Note that only results for the widest floating-point type on the system are given as 108narrower types have __zero_error. 109 110Note that errors grow as /a/ grows larger. 111 112Note also that the higher error rates for the 80 and 128 bit 113long double results are somewhat misleading: expected results that are 114zero at 64-bit double precision may be non-zero - but exceptionally small - 115with the larger exponent range of a long double. These results therefore 116reflect the more extreme nature of the tests conducted for these types. 117 118All values are in units of epsilon. 119 120[table_gamma_p] 121 122[table_gamma_q] 123 124[table_tgamma_lower] 125 126[table_tgamma_incomplete_] 127 128[h4 Testing] 129 130There are two sets of tests: spot tests compare values taken from 131[@http://functions.wolfram.com/GammaBetaErf/ Mathworld's online evaluator] 132with this implementation to perform a basic "sanity check". 133Accuracy tests use data generated at very high precision 134(using NTL's RR class set at 1000-bit precision) using this implementation 135with a very high precision 60-term __lanczos, and some but not all of the special 136case handling disabled. 137This is less than satisfactory: an independent method should really be used, 138but apparently a complete lack of such methods are available. We can't even use a deliberately 139naive implementation without special case handling since Legendre's continued fraction 140(see below) is unstable for small a and z. 141 142[h4 Implementation] 143 144These four functions share a common implementation since 145they are all related via: 146 1471) [equation igamma5] 148 1492) [equation igamma6] 150 1513) [equation igamma7] 152 153The lower incomplete gamma is computed from its series representation: 154 1554) [equation igamma8] 156 157Or by subtraction of the upper integral from either [Gamma](a) or 1 158when /x - (1/(3x)) > a and x > 1.1/. 159 160The upper integral is computed from Legendre's continued fraction representation: 161 1625) [equation igamma9] 163 164When /(x > 1.1)/ or by subtraction of the lower integral from either [Gamma](a) or 1 165when /x - (1/(3x)) < a/. 166 167For /x < 1.1/ computation of the upper integral is more complex as the continued 168fraction representation is unstable in this area. However there is another 169series representation for the lower integral: 170 1716) [equation igamma10] 172 173That lends itself to calculation of the upper integral via rearrangement 174to: 175 1767) [equation igamma11] 177 178Refer to the documentation for __powm1 and __tgamma1pm1 for details 179of their implementation. 180 181For /x < 1.1/ the crossover point where the result is ~0.5 no longer 182occurs for /x ~ y/. Using /x * 0.75 < a/ as the crossover criterion 183for /0.5 < x <= 1.1/ keeps the maximum value computed (whether 184it's the upper or lower interval) to around 0.75. Likewise for 185/x <= 0.5/ then using /-0.4 / log(x) < a/ as the crossover criterion 186keeps the maximum value computed to around 0.7 187(whether it's the upper or lower interval). 188 189There are two special cases used when a is an integer or half integer, 190and the crossover conditions listed above indicate that we should compute 191the upper integral Q. 192If a is an integer in the range /1 <= a < 30/ then the following 193finite sum is used: 194 1959) [equation igamma1f] 196 197While for half-integers in the range /0.5 <= a < 30/ then the 198following finite sum is used: 199 20010) [equation igamma2f] 201 202These are both more stable and more efficient than the continued fraction 203alternative. 204 205When the argument /a/ is large, and /x ~ a/ then the series (4) and continued 206fraction (5) above are very slow to converge. In this area an expansion due to 207Temme is used: 208 20911) [equation igamma16] 210 21112) [equation igamma17] 212 21313) [equation igamma18] 214 21514) [equation igamma19] 216 217The double sum is truncated to a fixed number of terms - to give a specific 218target precision - and evaluated as a polynomial-of-polynomials. There are 219versions for up to 128-bit long double precision: types requiring 220greater precision than that do not use these expansions. The 221coefficients C[sub k][super n] are computed in advance using the recurrence 222relations given by Temme. The zone where these expansions are used is 223 224 (a > 20) && (a < 200) && fabs(x-a)/a < 0.4 225 226And: 227 228 (a > 200) && (fabs(x-a)/a < 4.5/sqrt(a)) 229 230The latter range is valid for all types up to 128-bit long doubles, and 231is designed to ensure that the result is larger than 10[super -6], the 232first range is used only for types up to 80-bit long doubles. These 233domains are narrower than the ones recommended by either Temme or Didonato 234and Morris. However, using a wider range results in large and inexact 235(i.e. computed) values being passed to the `exp` and `erfc` functions 236resulting in significantly larger error rates. In other words there is a 237fine trade off here between efficiency and error. The current limits should 238keep the number of terms required by (4) and (5) to no more than ~20 239at double precision. 240 241For the normalised incomplete gamma functions, calculation of the 242leading power terms is central to the accuracy of the function. 243For smallish a and x combining 244the power terms with the __lanczos gives the greatest accuracy: 245 24615) [equation igamma12] 247 248In the event that this causes underflow/overflow then the exponent can 249be reduced by a factor of /a/ and brought inside the power term. 250 251When a and x are large, we end up with a very large exponent with a base 252near one: this will not be computed accurately via the pow function, 253and taking logs simply leads to cancellation errors. The worst of the 254errors can be avoided by using: 255 25616) [equation igamma13] 257 258when /a-x/ is small and a and x are large. There is still a subtraction 259and therefore some cancellation errors - but the terms are small so the absolute 260error will be small - and it is absolute rather than relative error that 261counts in the argument to the /exp/ function. Note that for sufficiently 262large a and x the errors will still get you eventually, although this does 263delay the inevitable much longer than other methods. Use of /log(1+x)-x/ here 264is inspired by Temme (see references below). 265 266[h4 References] 267 268* N. M. Temme, A Set of Algorithms for the Incomplete Gamma Functions, 269Probability in the Engineering and Informational Sciences, 8, 1994. 270* N. M. Temme, The Asymptotic Expansion of the Incomplete Gamma Functions, 271Siam J. Math Anal. Vol 10 No 4, July 1979, p757. 272* A. R. Didonato and A. H. Morris, Computation of the Incomplete Gamma 273Function Ratios and their Inverse. ACM TOMS, Vol 12, No 4, Dec 1986, p377. 274* W. Gautschi, The Incomplete Gamma Functions Since Tricomi, In Tricomi's Ideas 275and Contemporary Applied Mathematics, Atti dei Convegni Lincei, n. 147, 276Accademia Nazionale dei Lincei, Roma, 1998, pp. 203--237. 277[@http://citeseer.ist.psu.edu/gautschi98incomplete.html http://citeseer.ist.psu.edu/gautschi98incomplete.html] 278 279[endsect] [/section:igamma The Incomplete Gamma Function] 280 281[/ 282 Copyright 2006 John Maddock and Paul A. Bristow. 283 Distributed under the Boost Software License, Version 1.0. 284 (See accompanying file LICENSE_1_0.txt or copy at 285 http://www.boost.org/LICENSE_1_0.txt). 286] 287