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 27The [@http://en.wikipedia.org/wiki/Gamma_function lgamma function] is defined by: 28 29[equation lgamm1] 30 31The second form of the function takes a pointer to an integer, 32which if non-null is set on output to the sign of tgamma(z). 33 34[optional_policy] 35 36[graph lgamma] 37 38The return type of these functions is computed using the __arg_promotion_rules: 39the result is of type `double` if T is an integer type, or type T otherwise. 40 41[h4 Accuracy] 42 43The following table shows the peak errors (in units of epsilon) 44found on various platforms 45with various floating point types, along with comparisons to 46various other libraries. Unless otherwise specified any 47floating point type that is narrower than the one shown will have 48__zero_error. 49 50Note that while the relative errors near the positive roots of lgamma 51are very low, the lgamma function has an infinite number of irrational 52roots for negative arguments: very close to these negative roots only 53a low absolute error can be guaranteed. 54 55[table_lgamma] 56 57The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision, 58and GCC-7.1/Ubuntu for `long double` and `__float128`. 59 60[graph lgamma__double] 61 62[graph lgamma__80_bit_long_double] 63 64[graph lgamma____float128] 65 66[h4 Testing] 67 68The main tests for this function involve comparisons against the logs of 69the factorials which can be independently calculated to very high accuracy. 70 71Random tests in key problem areas are also used. 72 73[h4 Implementation] 74 75The generic version of this function is implemented using Sterling's approximation 76for large arguments: 77 78[equation gamma6] 79 80For small arguments, the logarithm of tgamma is used. 81 82For negative /z/ the logarithm version of the 83reflection formula is used: 84 85[equation lgamm3] 86 87For types of known precision, the __lanczos is used, a traits class 88`boost::math::lanczos::lanczos_traits` maps type T to an appropriate 89approximation. The logarithmic version of the __lanczos is: 90 91[equation lgamm4] 92 93Where L[sub e,g] is the Lanczos sum, scaled by e[super g]. 94 95As before the reflection formula is used for /z < 0/. 96 97When z is very near 1 or 2, then the logarithmic version of the __lanczos 98suffers very badly from cancellation error: indeed for values sufficiently 99close to 1 or 2, arbitrarily large relative errors can be obtained (even though 100the absolute error is tiny). 101 102For types with up to 113 bits of precision 103(up to and including 128-bit long doubles), root-preserving 104rational approximations [jm_rationals] are used 105over the intervals [1,2] and [2,3]. Over the interval [2,3] the approximation 106form used is: 107 108 lgamma(z) = (z-2)(z+1)(Y + R(z-2)); 109 110Where Y is a constant, and R(z-2) is the rational approximation: optimised 111so that its absolute error is tiny compared to Y. In addition, small values of z greater 112than 3 can handled by argument reduction using the recurrence relation: 113 114 lgamma(z+1) = log(z) + lgamma(z); 115 116Over 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 120Once again Y is a constant, and R(z-1) is optimised for low absolute error 121compared to Y. For z > 1.5 the above form wouldn't converge to a 122minimax solution but this similar form does: 123 124 lgamma(z) = (2-z)(1-z)(Y + R(2-z)); 125 126Finally for z < 1 the recurrence relation can be used to move to z > 1: 127 128 lgamma(z) = lgamma(z+1) - log(z); 129 130Note that while this involves a subtraction, it appears not 131to suffer from cancellation error: as z decreases from 1 132the `-log(z)` term grows positive much more 133rapidly than the `lgamma(z+1)` term becomes negative. So in this 134specific case, significant digits are preserved, rather than cancelled. 135 136For other types which do have a __lanczos defined for them 137the current solution is as follows: imagine we 138balance the two terms in the __lanczos by dividing the power term by its value 139at /z = 1/, and then multiplying the Lanczos coefficients by the same value. 140Now each term will take the value 1 at /z = 1/ and we can rearrange the power terms 141in 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 143a new summation that can be also be fed into log1p. Crucially, all of the 144terms tend to zero, as /z -> 1/: 145 146[equation lgamm5] 147 148The C[sub k] terms in the above are the same as in the __lanczos. 149 150A 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] 162