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_pomotion_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 the
  54 __gsl, __glibc, __hpc and
  55 __cephes libraries. Unless otherwise specified any
  56 floating point type that is narrower than the one shown will have
  57 __zero_error.
  58
  59 Note that while the relative errors near the positive roots of lgamma
  60 are very low, the lgamma function has an infinite number of irrational
  61 roots for negative arguments: very close to these negative roots only
  62 a low absolute error can be guaranteed.
  63
  64 [table
  65 [[Significand Size] [Platform and Compiler] [Factorials and Half factorials] [Values Near Zero] [Values Near 1 or 2] [Values Near a Negative Pole]]
  66 [[53]            [Win32 Visual C++ 8]
  67 [Peak=0.88 Mean=0.14
  68
  69 (GSL=33) (__cephes=1.5)]
  70 [Peak=0.96 Mean=0.46
  71
  72 (GSL=5.2) (__cephes=1.1)]
  73 [Peak=0.86 Mean=0.46
  74
  75 (GSL=1168) (__cephes~500000)]
  76 [Peak=4.2 Mean=1.3
  77
  78 (GSL=25) (__cephes=1.6)] ]
  79 [[64]            [Linux IA32 / GCC]
  80 [Peak=1.9 Mean=0.43
  81
  82 (__glibc Peak=1.7 Mean=0.49)]
  83 [Peak=1.4 Mean=0.57
  84
  85 (__glibc Peak= 0.96 Mean=0.54)]
  86 [Peak=0.86 Mean=0.35
  87
  88 (__glibc Peak=0.74 Mean=0.26)]
  89
  90 [Peak=6.0 Mean=1.8
  91
  92 (__glibc Peak=3.0 Mean=0.86)]           ]
  93 [[64]            [Linux IA64 / GCC]
  94 [Peak=0.99 Mean=0.12
  95
  96 (__glibc Peak 0)]
  97
  98 [Pek=1.2 Mean=0.6
  99
 100 (__glibc Peak 0)]
 101 [Peak=0.86 Mean=0.16
 102
 103 (__glibc Peak 0)]
 104 [Peak=2.3 Mean=0.69
 105
 106 (__glibc Peak 0)]           ]
 107 [[113]           [HPUX IA64, aCC A.06.06]
 108 [Peak=0.96 Mean=0.13
 109
 110 (__hpc Peak 0)]
 111 [Peak=0.99 Mean=0.53
 112
 113 (__hpc Peak 0)]
 114 [Peak=0.9 Mean=0.4
 115
 116 (__hpc Peak 0)]
 117 [Peak=3.0 Mean=0.9
 118
 119 (__hpc Peak 0)]  ]
 120 ]
 121
 122 [h4 Testing]
 123
 124 The main tests for this function involve comparisons against the logs of
 125 the factorials which can be independently calculated to very high accuracy.
 126
 127 Random tests in key problem areas are also used.
 128
 129 [h4 Implementation]
 130
 131 The generic version of this function is implemented using Sterling's approximation
 132 for large arguments:
 133
 134 [equation gamma6]
 135
 136 For small arguments, the logarithm of tgamma is used.
 137
 138 For negative /z/ the logarithm version of the
 139 reflection formula is used:
 140
 141 [equation lgamm3]
 142
 143 For types of known precision, the __lanczos is used, a traits class
 144 `boost::math::lanczos::lanczos_traits` maps type T to an appropriate
 145 approximation.  The logarithmic version of the __lanczos is:
 146
 147 [equation lgamm4]
 148
 149 Where L[sub e,g][space] is the Lanczos sum, scaled by e[super g].
 150
 151 As before the reflection formula is used for /z < 0/.
 152
 153 When z is very near 1 or 2, then the logarithmic version of the __lanczos
 154 suffers very badly from cancellation error: indeed for values sufficiently
 155 close to 1 or 2, arbitrarily large relative errors can be obtained (even though
 156 the absolute error is tiny).
 157
 158 For types with up to 113 bits of precision
 159 (up to and including 128-bit long doubles), root-preserving
 160 rational approximations [jm_rationals] are used
 161 over the intervals [1,2] and [2,3].  Over the interval [2,3] the approximation
 162 form used is:
 163
 164    lgamma(z) = (z-2)(z+1)(Y + R(z-2));
 165
 166 Where Y is a constant, and R(z-2) is the rational approximation: optimised
 167 so that it's absolute error is tiny compared to Y.  In addition
 168 small values of z greater
 169 than 3 can handled by argument reduction using the recurrence relation:
 170
 171    lgamma(z+1) = log(z) + lgamma(z);
 172
 173 Over the interval [1,2] two approximations have to be used, one for small z uses:
 174
 175    lgamma(z) = (z-1)(z-2)(Y + R(z-1));
 176
 177 Once again Y is a constant, and R(z-1) is optimised for low absolute error
 178 compared to Y.  For z > 1.5 the above form wouldn't converge to a
 179 minimax solution but this similar form does:
 180
 181    lgamma(z) = (2-z)(1-z)(Y + R(2-z));
 182
 183 Finally for z < 1 the recurrence relation can be used to move to z > 1:
 184
 185    lgamma(z) = lgamma(z+1) - log(z);
 186
 187 Note that while this involves a subtraction, it appears not
 188 to suffer from cancellation error: as z decreases from 1
 189 the `-log(z)` term grows positive much more
 190 rapidly than the `lgamma(z+1)` term becomes negative.  So in this
 191 specific case, significant digits are preserved, rather than cancelled.
 192
 193 For other types which do have a __lanczos defined for them
 194 the current solution is as follows: imagine we
 195 balance the two terms in the __lanczos by dividing the power term by its value
 196 at /z = 1/, and then multiplying the Lanczos coefficients by the same value.
 197 Now each term will take the value 1 at /z = 1/ and we can rearrange the power terms
 198 in terms of log1p.  Likewise if we subtract 1 from the Lanczos sum part
 199 (algebraically, by subtracting the value of each term at /z = 1/), we obtain
 200 a new summation that can be also be fed into log1p.  Crucially, all of the
 201 terms tend to zero, as /z -> 1/:
 202
 203 [equation lgamm5]
 204
 205 The C[sub k][space] terms in the above are the same as in the __lanczos.
 206
 207 A similar rearrangement can be performed at /z = 2/:
 208
 209 [equation lgamm6]
 210
 211 [endsect][/section:lgamma The Log Gamma Function]
 212
 213 [/
 214   Copyright 2006 John Maddock and Paul A. Bristow.
 215   Distributed under the Boost Software License, Version 1.0.
 216   (See accompanying file LICENSE_1_0.txt or copy at
 217   http://www.boost.org/LICENSE_1_0.txt).
 218 ]