libs/math/doc/sf/igamma.qbk

   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
  39 There are four [@http://mathworld.wolfram.com/IncompleteGammaFunction.html
  40 incomplete gamma functions]:
  41 two are normalised versions (also known as /regularized/ incomplete gamma functions)
  42 that return values in the range [0, 1], and two are non-normalised and
  43 return values in the range [0, [Gamma](a)].  Users interested in statistical
  44 applications should use the
  45 [@http://mathworld.wolfram.com/RegularizedGammaFunction.html normalised versions (`gamma_p` and `gamma_q`)].
  46
  47 All of these functions require /a > 0/ and /z >= 0/, otherwise they return
  48 the result of __domain_error.
  49
  50 [optional_policy]
  51
  52 The return type of these functions is computed using the __arg_promotion_rules
  53 when 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
  61 Returns the normalised lower incomplete gamma function of a and z:
  62
  63 [equation igamma4]
  64
  65 This 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
  75 Returns the normalised upper incomplete gamma function of a and z:
  76
  77 [equation igamma3]
  78
  79 This 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
  89 Returns 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
  99 Returns the full (non-normalised) upper incomplete gamma function of a and z:
 100
 101 [equation igamma1]
 102
 103 [h4 Accuracy]
 104
 105 The following tables give peak and mean relative errors in over various domains of
 106 a and z, along with comparisons to the __gsl and __cephes libraries.
 107 Note that only results for the widest floating-point type on the system are given as
 108 narrower types have __zero_error.
 109
 110 Note that errors grow as /a/ grows larger.
 111
 112 Note also that the higher error rates for the 80 and 128 bit
 113 long double results are somewhat misleading: expected results that are
 114 zero at 64-bit double precision may be non-zero - but exceptionally small -
 115 with the larger exponent range of a long double.  These results therefore
 116 reflect the more extreme nature of the tests conducted for these types.
 117
 118 All 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
 130 There are two sets of tests: spot tests compare values taken from
 131 [@http://functions.wolfram.com/GammaBetaErf/ Mathworld's online evaluator]
 132 with this implementation to perform a basic "sanity check".
 133 Accuracy tests use data generated at very high precision
 134 (using NTL's RR class set at 1000-bit precision) using this implementation
 135 with a very high precision 60-term __lanczos, and some but not all of the special
 136 case handling disabled.
 137 This is less than satisfactory: an independent method should really be used,
 138 but apparently a complete lack of such methods are available.  We can't even use a deliberately
 139 naive 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
 144 These four functions share a common implementation since
 145 they are all related via:
 146
 147 1) [equation igamma5]
 148
 149 2) [equation igamma6]
 150
 151 3) [equation igamma7]
 152
 153 The lower incomplete gamma is computed from its series representation:
 154
 155 4) [equation igamma8]
 156
 157 Or by subtraction of the upper integral from either [Gamma](a) or 1
 158 when /x - (1/(3x)) > a and x > 1.1/.
 159
 160 The upper integral is computed from Legendre's continued fraction representation:
 161
 162 5) [equation igamma9]
 163
 164 When /(x > 1.1)/ or by subtraction of the lower integral from either [Gamma](a) or 1
 165 when /x - (1/(3x))  < a/.
 166
 167 For /x < 1.1/ computation of the upper integral is more complex as the continued
 168 fraction representation is unstable in this area.  However there is another
 169 series representation for the lower integral:
 170
 171 6) [equation igamma10]
 172
 173 That lends itself to calculation of the upper integral via rearrangement
 174 to:
 175
 176 7) [equation igamma11]
 177
 178 Refer to the documentation for __powm1 and __tgamma1pm1 for details
 179 of their implementation.
 180
 181 For /x < 1.1/ the crossover point where the result is ~0.5 no longer
 182 occurs for /x ~ y/.  Using /x * 0.75 < a/ as the crossover criterion
 183 for /0.5 < x <= 1.1/ keeps the maximum value computed (whether
 184 it'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
 186 keeps the maximum value computed to around 0.7
 187 (whether it's the upper or lower interval).
 188
 189 There are two special cases used when a is an integer or half integer,
 190 and the crossover conditions listed above indicate that we should compute
 191 the upper integral Q.
 192 If a is an integer in the range /1 <= a < 30/ then the following
 193 finite sum is used:
 194
 195 9) [equation igamma1f]
 196
 197 While for half-integers in the range /0.5 <= a < 30/ then the
 198 following finite sum is used:
 199
 200 10) [equation igamma2f]
 201
 202 These are both more stable and more efficient than the continued fraction
 203 alternative.
 204
 205 When the argument /a/ is large, and /x ~ a/ then the series (4) and continued
 206 fraction (5) above are very slow to converge.  In this area an expansion due to
 207 Temme is used:
 208
 209 11) [equation igamma16]
 210
 211 12) [equation igamma17]
 212
 213 13) [equation igamma18]
 214
 215 14) [equation igamma19]
 216
 217 The double sum is truncated to a fixed number of terms - to give a specific
 218 target precision - and evaluated as a polynomial-of-polynomials.  There are
 219 versions for up to 128-bit long double precision: types requiring
 220 greater precision than that do not use these expansions.  The
 221 coefficients C[sub k][super n] are computed in advance using the recurrence
 222 relations given by Temme.  The zone where these expansions are used is
 223
 224    (a > 20) && (a < 200) && fabs(x-a)/a < 0.4
 225
 226 And:
 227
 228    (a > 200) && (fabs(x-a)/a < 4.5/sqrt(a))
 229
 230 The latter range is valid for all types up to 128-bit long doubles, and
 231 is designed to ensure that the result is larger than 10[super -6], the
 232 first range is used only for types up to 80-bit long doubles.  These
 233 domains are narrower than the ones recommended by either Temme or Didonato
 234 and Morris.  However, using a wider range results in large and inexact
 235 (i.e. computed) values being passed to the `exp` and `erfc` functions
 236 resulting in significantly larger error rates.  In other words there is a
 237 fine trade off here between efficiency and error.  The current limits should
 238 keep the number of terms required by (4) and (5) to no more than ~20
 239 at double precision.
 240
 241 For the normalised incomplete gamma functions, calculation of the
 242 leading power terms is central to the accuracy of the function.
 243 For smallish a and x combining
 244 the power terms with the __lanczos gives the greatest accuracy:
 245
 246 15) [equation igamma12]
 247
 248 In the event that this causes underflow/overflow then the exponent can
 249 be reduced by a factor of /a/ and brought inside the power term.
 250
 251 When a and x are large, we end up with a very large exponent with a base
 252 near one: this will not be computed accurately via the pow function,
 253 and taking logs simply leads to cancellation errors.  The worst of the
 254 errors can be avoided by using:
 255
 256 16) [equation igamma13]
 257
 258 when /a-x/ is small and a and x are large.  There is still a subtraction
 259 and therefore some cancellation errors - but the terms are small so the absolute
 260 error will be small - and it is absolute rather than relative error that
 261 counts in the argument to the /exp/ function.  Note that for sufficiently
 262 large a and x the errors will still get you eventually, although this does
 263 delay the inevitable much longer than other methods.  Use of /log(1+x)-x/ here
 264 is inspired by Temme (see references below).
 265
 266 [h4 References]
 267
 268 * N. M. Temme, A Set of Algorithms for the Incomplete Gamma Functions,
 269 Probability in the Engineering and Informational Sciences, 8, 1994.
 270 * N. M. Temme, The Asymptotic Expansion of the Incomplete Gamma Functions,
 271 Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
 272 * A. R. Didonato and A. H. Morris, Computation of the Incomplete Gamma
 273 Function 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
 275 and Contemporary Applied Mathematics, Atti dei Convegni Lincei, n. 147,
 276 Accademia 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 ]