]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ] |