ceph/src/boost/libs/math/doc/sf/lgamma.qbk

   1 [section:lgamma Log Gamma]
   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 T>
  12    ``__sf_result`` lgamma(T z);
  13
  14    template <class T, class ``__Policy``>
  15    ``__sf_result`` lgamma(T z, const ``__Policy``&);
  16
  17    template <class T>
  18    ``__sf_result`` lgamma(T z, int* sign);
  19
  20    template <class T, class ``__Policy``>
  21    ``__sf_result`` lgamma(T z, int* sign, const ``__Policy``&);
  22
  23    }} // namespaces
  24
  25 [h4 Description]
  26
  27 The [@http://en.wikipedia.org/wiki/Gamma_function lgamma function] is defined by:
  28
  29 [equation lgamm1]
  30
  31 The second form of the function takes a pointer to an integer,
  32 which if non-null is set on output to the sign of tgamma(z).
  33
  34 [optional_policy]
  35
  36 [graph lgamma]
  37
  38 There are effectively two versions of this function internally: a fully
  39 generic version that is slow, but reasonably accurate, and a much more
  40 efficient approximation that is used where the number of digits in the significand
  41 of T correspond to a certain __lanczos.  In practice, any built-in
  42 floating-point type you will encounter has an appropriate __lanczos
  43 defined for it.  It is also possible, given enough machine time, to generate
  44 further __lanczos's using the program libs/math/tools/lanczos_generator.cpp.
  45
  46 The return type of these functions is computed using the __arg_promotion_rules:
  47 the result is of type `double` if T is an integer type, or type T otherwise.
  48
  49 [h4 Accuracy]
  50
  51 The following table shows the peak errors (in units of epsilon)
  52 found on various platforms
  53 with various floating point types, along with comparisons to
  54 various other libraries. Unless otherwise specified any
  55 floating point type that is narrower than the one shown will have
  56 __zero_error.
  57
  58 Note that while the relative errors near the positive roots of lgamma
  59 are very low, the lgamma function has an infinite number of irrational
  60 roots for negative arguments: very close to these negative roots only
  61 a low absolute error can be guaranteed.
  62
  63 [table_lgamma]
  64
  65 [h4 Testing]
  66
  67 The main tests for this function involve comparisons against the logs of
  68 the factorials which can be independently calculated to very high accuracy.
  69
  70 Random tests in key problem areas are also used.
  71
  72 [h4 Implementation]
  73
  74 The generic version of this function is implemented using Sterling's approximation
  75 for large arguments:
  76
  77 [equation gamma6]
  78
  79 For small arguments, the logarithm of tgamma is used.
  80
  81 For negative /z/ the logarithm version of the
  82 reflection formula is used:
  83
  84 [equation lgamm3]
  85
  86 For types of known precision, the __lanczos is used, a traits class
  87 `boost::math::lanczos::lanczos_traits` maps type T to an appropriate
  88 approximation.  The logarithmic version of the __lanczos is:
  89
  90 [equation lgamm4]
  91
  92 Where L[sub e,g][space] is the Lanczos sum, scaled by e[super g].
  93
  94 As before the reflection formula is used for /z < 0/.
  95
  96 When z is very near 1 or 2, then the logarithmic version of the __lanczos
  97 suffers very badly from cancellation error: indeed for values sufficiently
  98 close to 1 or 2, arbitrarily large relative errors can be obtained (even though
  99 the absolute error is tiny).
 100
 101 For types with up to 113 bits of precision
 102 (up to and including 128-bit long doubles), root-preserving
 103 rational approximations [jm_rationals] are used
 104 over the intervals [1,2] and [2,3].  Over the interval [2,3] the approximation
 105 form used is:
 106
 107    lgamma(z) = (z-2)(z+1)(Y + R(z-2));
 108
 109 Where Y is a constant, and R(z-2) is the rational approximation: optimised
 110 so that it's absolute error is tiny compared to Y.  In addition
 111 small values of z greater
 112 than 3 can handled by argument reduction using the recurrence relation:
 113
 114    lgamma(z+1) = log(z) + lgamma(z);
 115
 116 Over the interval [1,2] two approximations have to be used, one for small z uses:
 117
 118    lgamma(z) = (z-1)(z-2)(Y + R(z-1));
 119
 120 Once again Y is a constant, and R(z-1) is optimised for low absolute error
 121 compared to Y.  For z > 1.5 the above form wouldn't converge to a
 122 minimax solution but this similar form does:
 123
 124    lgamma(z) = (2-z)(1-z)(Y + R(2-z));
 125
 126 Finally for z < 1 the recurrence relation can be used to move to z > 1:
 127
 128    lgamma(z) = lgamma(z+1) - log(z);
 129
 130 Note that while this involves a subtraction, it appears not
 131 to suffer from cancellation error: as z decreases from 1
 132 the `-log(z)` term grows positive much more
 133 rapidly than the `lgamma(z+1)` term becomes negative.  So in this
 134 specific case, significant digits are preserved, rather than cancelled.
 135
 136 For other types which do have a __lanczos defined for them
 137 the current solution is as follows: imagine we
 138 balance the two terms in the __lanczos by dividing the power term by its value
 139 at /z = 1/, and then multiplying the Lanczos coefficients by the same value.
 140 Now each term will take the value 1 at /z = 1/ and we can rearrange the power terms
 141 in terms of log1p.  Likewise if we subtract 1 from the Lanczos sum part
 142 (algebraically, by subtracting the value of each term at /z = 1/), we obtain
 143 a new summation that can be also be fed into log1p.  Crucially, all of the
 144 terms tend to zero, as /z -> 1/:
 145
 146 [equation lgamm5]
 147
 148 The C[sub k][space] terms in the above are the same as in the __lanczos.
 149
 150 A similar rearrangement can be performed at /z = 2/:
 151
 152 [equation lgamm6]
 153
 154 [endsect][/section:lgamma The Log Gamma Function]
 155
 156 [/
 157   Copyright 2006 John Maddock and Paul A. Bristow.
 158   Distributed under the Boost Software License, Version 1.0.
 159   (See accompanying file LICENSE_1_0.txt or copy at
 160   http://www.boost.org/LICENSE_1_0.txt).
 161 ]