igamma.qbk

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

[section:igamma Incomplete Gamma Functions]

[h4 Synopsis]

``
#include <boost/math/special_functions/gamma.hpp>
``

   namespace boost{ namespace math{
   
   template <class T1, class T2>
   ``__sf_result`` gamma_p(T1 a, T2 z);
   
   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` gamma_p(T1 a, T2 z, const ``__Policy``&);
   
   template <class T1, class T2>
   ``__sf_result`` gamma_q(T1 a, T2 z);
   
   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` gamma_q(T1 a, T2 z, const ``__Policy``&);
   
   template <class T1, class T2>
   ``__sf_result`` tgamma_lower(T1 a, T2 z);
   
   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` tgamma_lower(T1 a, T2 z, const ``__Policy``&);
   
   template <class T1, class T2>
   ``__sf_result`` tgamma(T1 a, T2 z);
   
   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` tgamma(T1 a, T2 z, const ``__Policy``&);
   
   }} // namespaces
   
[h4 Description]

There are four [@http://mathworld.wolfram.com/IncompleteGammaFunction.html 
incomplete gamma functions]:
two are normalised versions (also known as /regularized/ incomplete gamma functions)
that return values in the range [0, 1], and two are non-normalised and
return values in the range [0, [Gamma](a)].  Users interested in statistical
applications should use the
[@http://mathworld.wolfram.com/RegularizedGammaFunction.html normalised versions (gamma_p and gamma_q)].

All of these functions require /a > 0/ and /z >= 0/, otherwise they return
the result of __domain_error.

[optional_policy]

The return type of these functions is computed using the __arg_pomotion_rules
when T1 and T2 are different types, otherwise the return type is simply T1.

   template <class T1, class T2>
   ``__sf_result`` gamma_p(T1 a, T2 z);
   
   template <class T1, class T2, class Policy>
   ``__sf_result`` gamma_p(T1 a, T2 z, const ``__Policy``&);
   
Returns the normalised lower incomplete gamma function of a and z:

[equation igamma4]

This function changes rapidly from 0 to 1 around the point z == a:

[graph gamma_p]

   template <class T1, class T2>
   ``__sf_result`` gamma_q(T1 a, T2 z);

   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` gamma_q(T1 a, T2 z, const ``__Policy``&);

Returns the normalised upper incomplete gamma function of a and z:

[equation igamma3]

This function changes rapidly from 1 to 0 around the point z == a:

[graph gamma_q]

   template <class T1, class T2>
   ``__sf_result`` tgamma_lower(T1 a, T2 z);

   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` tgamma_lower(T1 a, T2 z, const ``__Policy``&);

Returns the full (non-normalised) lower incomplete gamma function of a and z:

[equation igamma2]

   template <class T1, class T2>
   ``__sf_result`` tgamma(T1 a, T2 z);

   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` tgamma(T1 a, T2 z, const ``__Policy``&);

Returns the full (non-normalised) upper incomplete gamma function of a and z:

[equation igamma1]

[h4 Accuracy]

The following tables give peak and mean relative errors in over various domains of
a and z, along with comparisons to the __gsl and __cephes libraries.
Note that only results for the widest floating point type on the system are given as
narrower types have __zero_error.

Note that errors grow as /a/ grows larger.

Note also that the higher error rates for the 80 and 128 bit 
long double results are somewhat misleading: expected results that are 
zero at 64-bit double precision may be non-zero - but exceptionally small -
with the larger exponent range of a long double.  These results therefore
reflect the more extreme nature of the tests conducted for these types.

All values are in units of epsilon.

[table Errors In the Function gamma_p(a,z)
[[Significand Size] [Platform and Compiler] 
   [0.5 < a < 100

   and

   0.01*a < z < 100*a]  
      [1x10[super -12] < a < 5x10[super -2]

      and

      0.01*a < z < 100*a] 
         [1e-6 < a < 1.7x10[super 6]

         and

         1 < z < 100*a]]
[[53] [Win32, Visual C++ 8] 
   [Peak=36 Mean=9.1

   (GSL Peak=342 Mean=46)

   (__cephes Peak=491 Mean=102)]  
      [Peak=4.5 Mean=1.4

      (GSL Peak=4.8 Mean=0.76)

      (__cephes Peak=21 Mean=5.6)] 
         [Peak=244 Mean=21

         (GSL Peak=1022 Mean=1054)

         (__cephes Peak~8x10[super 6] Mean~7x10[super 4])]]
[[64] [RedHat Linux IA32, gcc-3.3] [Peak=241 Mean=36]  [Peak=4.7 Mean=1.5] [Peak~30,220 Mean=1929]]
[[64] [Redhat Linux IA64, gcc-3.4] [Peak=41 Mean=10]  [Peak=4.7 Mean=1.4] [Peak~30,790 Mean=1864]]
[[113] [HPUX IA64, aCC A.06.06] [Peak=40.2 Mean=10.2]  [Peak=5 Mean=1.6] [Peak=5,476 Mean=440]]
]

[table Errors In the Function gamma_q(a,z)
[[Significand Size] [Platform and Compiler] 
[0.5 < a < 100

and

0.01*a < z < 100*a]  
[1x10[super -12] < a < 5x10[super -2]

and

0.01*a < z < 100*a] [1x10[super -6] < a < 1.7x10[super 6]

and

1 < z < 100*a]]
[[53] [Win32, Visual C++ 8] [Peak=28.3 Mean=7.2

(GSL Peak=201 Mean=13)

(__cephes Peak=556 Mean=97)]  [Peak=4.8 Mean=1.6

(GSL Peak~1.3x10[super 10] Mean=1x10[super +9])

(__cephes Peak~3x10[super 11] Mean=4x10[super 10])] [Peak=469 Mean=33

(GSL Peak=27,050 Mean=2159)

(__cephes Peak~8x10[super 6] Mean~7x10[super 5])]]
[[64] [RedHat Linux IA32, gcc-3.3] [Peak=280 Mean=33]  [Peak=4.1 Mean=1.6] [Peak=11,490 Mean=732]]
[[64] [Redhat Linux IA64, gcc-3.4] [Peak=32 Mean=9.4]  [Peak=4.7 Mean=1.5] [Peak=6815 Mean=414]]
[[113] [HPUX IA64, aCC A.06.06] [Peak=37 Mean=10]  [Peak=11.2 Mean=2.0] [Peak=4,999 Mean=298]]
]

[table Errors In the Function tgamma_lower(a,z)
[[Significand Size] [Platform and Compiler] [0.5 < a < 100

and

0.01*a < z < 100*a]  [1x10[super -12] < a < 5x10[super -2]

and

0.01*a < z < 100*a]]
[[53] [Win32, Visual C++ 8] [Peak=5.5 Mean=1.4]  [Peak=3.6 Mean=0.78]]
[[64] [RedHat Linux IA32, gcc-3.3] [Peak=402 Mean=79]  [Peak=3.4 Mean=0.8]]
[[64] [Redhat Linux IA64, gcc-3.4] [Peak=6.8 Mean=1.4]  [Peak=3.4 Mean=0.78]]
[[113] [HPUX IA64, aCC A.06.06] [Peak=6.1 Mean=1.8]  [Peak=3.7 Mean=0.89]]
]

[table Errors In the Function tgamma(a,z)
[[Significand Size] [Platform and Compiler] [0.5 < a < 100

and

0.01*a < z < 100*a]  [1x10[super -12] < a < 5x10[super -2]

and

0.01*a < z < 100*a]]
[[53] [Win32, Visual C++ 8] [Peak=5.9 Mean=1.5]  [Peak=1.8 Mean=0.6]]
[[64] [RedHat Linux IA32, gcc-3.3] [Peak=596 Mean=116]  [Peak=3.2 Mean=0.84]]
[[64] [Redhat Linux IA64, gcc-3.4.4] [Peak=40.2 Mean=2.5]  [Peak=3.2 Mean=0.8]]
[[113] [HPUX IA64, aCC A.06.06] [Peak=364 Mean=17.6]  [Peak=12.7 Mean=1.8]]
]

[h4 Testing]

There are two sets of tests: spot tests compare values taken from
[@http://functions.wolfram.com/GammaBetaErf/ Mathworld's online evaluator]
with this implementation to perform a basic "sanity check".
Accuracy tests use data generated at very high precision
(using NTL's RR class set at 1000-bit precision) using this implementation 
with a very high precision 60-term __lanczos, and some but not all of the special
case handling disabled.
This is less than satisfactory: an independent method should really be used,
but apparently a complete lack of such methods are available.  We can't even use a deliberately
naive implementation without special case handling since Legendre's continued fraction
(see below) is unstable for small a and z.

[h4 Implementation]

These four functions share a common implementation since
they are all related via:

1) [equation igamma5]

2) [equation igamma6]

3) [equation igamma7]

The lower incomplete gamma is computed from its series representation:

4) [equation igamma8]

Or by subtraction of the upper integral from either [Gamma](a) or 1
when /x > a and x > 1.1/.

The upper integral is computed from Legendre's continued fraction representation:

5) [equation igamma9]

When /x > 1.1/ or by subtraction of the lower integral from either [Gamma](a) or 1
when /x < a/.

For /x < 1.1/ computation of the upper integral is more complex as the continued 
fraction representation is unstable in this area.  However there is another 
series representation for the lower integral:

6) [equation igamma10]

That lends itself to calculation of the upper integral via rearrangement
to:

7) [equation igamma11]

Refer to the documentation for __powm1 and __tgamma1pm1 for details
of their implementation.  Note however that the precision of __tgamma1pm1
is capped to either around 35 digits, or to that of the __lanczos associated with
type T - if there is one - whichever of the two is the greater.  
That therefore imposes a similar limit on the precision of this 
function in this region.

For /x < 1.1/ the crossover point where the result is ~0.5 no longer
occurs for /x ~ y/.  Using /x * 1.1 < a/ as the crossover criterion
for /0.5 < x <= 1.1/ keeps the maximum value computed (whether
it's the upper or lower interval) to around 0.6.   Likewise for
/x <= 0.5/ then using /-0.4 / log(x) < a/ as the crossover criterion
keeps the maximum value computed to around 0.7
(whether it's the upper or lower interval).

There are two special cases used when a is an integer or half integer,
and the crossover conditions listed above indicate that we should compute
the upper integral Q.
If a is an integer in the range /1 <= a < 30/ then the following 
finite sum is used:

9) [equation igamma1f]

While for half integers in the range /0.5 <= a < 30/ then the
following finite sum is used:

10) [equation igamma2f]

These are both more stable and more efficient than the continued fraction
alternative.

When the argument /a/ is large, and /x ~ a/ then the series (4) and continued 
fraction (5) above are very slow to converge.  In this area an expansion due to
Temme is used:

11) [equation igamma16]

12) [equation igamma17]

13) [equation igamma18]

14) [equation igamma19]

The double sum is truncated to a fixed number of terms - to give a specific
target precision - and evaluated as a polynomial-of-polynomials.  There are 
versions for up to 128-bit long double precision: types requiring
greater precision than that do not use these expansions.  The
coefficients C[sub k][super n] are computed in advance using the recurrence
relations given by Temme.  The zone where these expansions are used is

   (a > 20) && (a < 200) && fabs(x-a)/a < 0.4
   
And:

   (a > 200) && (fabs(x-a)/a < 4.5/sqrt(a))
   
The latter range is valid for all types up to 128-bit long doubles, and
is designed to ensure that the result is larger than 10[super -6], the 
first range is used only for types up to 80-bit long doubles.  These
domains are narrower than the ones recommended by either Temme or Didonato
and Morris.  However, using a wider range results in large and inexact
(i.e. computed) values being passed to the `exp` and `erfc` functions
resulting in significantly larger error rates.  In other words there is a
fine trade off here between efficiency and error.  The current limits should
keep the number of terms required by (4) and (5) to no more than ~20
at double precision.

For the normalised incomplete gamma functions, calculation of the 
leading power terms is central to the accuracy of the function.
For smallish a and x combining
the power terms with the __lanczos gives the greatest accuracy:

15) [equation igamma12]

In the event that this causes underflow/overflow then the exponent can 
be reduced by a factor of /a/ and brought inside the power term.

When a and x are large, we end up with a very large exponent with a base
near one: this will not be computed accurately via the pow function,
and taking logs simply leads to cancellation errors.  The worst of the
errors can be avoided by using:

16) [equation igamma13]

when /a-x/ is small and a and x are large.  There is still a subtraction
and therefore some cancellation errors - but the terms are small so the absolute
error will be small - and it is absolute rather than relative error that 
counts in the argument to the /exp/ function.  Note that for sufficiently
large a and x the errors will still get you eventually, although this does
delay the inevitable much longer than other methods.  Use of /log(1+x)-x/ here
is inspired by Temme (see references below).

[h4 References]

* N. M. Temme, A Set of Algorithms for the Incomplete Gamma Functions,
Probability in the Engineering and Informational Sciences, 8, 1994.
* N. M. Temme, The Asymptotic Expansion of the Incomplete Gamma Functions,
Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
* A. R. Didonato and A. H. Morris, Computation of the Incomplete Gamma 
Function Ratios and their Inverse.  ACM TOMS, Vol 12, No 4, Dec 1986, p377.
* W. Gautschi, The Incomplete Gamma Functions Since Tricomi, In Tricomi's Ideas 
and Contemporary Applied Mathematics, Atti dei Convegni Lincei, n. 147, 
Accademia Nazionale dei Lincei, Roma, 1998, pp. 203--237. 
[@http://citeseer.ist.psu.edu/gautschi98incomplete.html http://citeseer.ist.psu.edu/gautschi98incomplete.html]

[endsect][/section:igamma The Incomplete Gamma Function]

[/ 
  Copyright 2006 John Maddock and Paul A. Bristow.
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or copy at
  http://www.boost.org/LICENSE_1_0.txt).
]