]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | [section:tgamma 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`` tgamma(T z); | |
13 | ||
14 | template <class T, class ``__Policy``> | |
15 | ``__sf_result`` tgamma(T z, const ``__Policy``&); | |
16 | ||
17 | template <class T> | |
18 | ``__sf_result`` tgamma1pm1(T dz); | |
19 | ||
20 | template <class T, class ``__Policy``> | |
21 | ``__sf_result`` tgamma1pm1(T dz, const ``__Policy``&); | |
22 | ||
23 | }} // namespaces | |
24 | ||
25 | [h4 Description] | |
26 | ||
27 | template <class T> | |
28 | ``__sf_result`` tgamma(T z); | |
29 | ||
30 | template <class T, class ``__Policy``> | |
31 | ``__sf_result`` tgamma(T z, const ``__Policy``&); | |
32 | ||
33 | Returns the "true gamma" (hence name tgamma) of value z: | |
34 | ||
35 | [equation gamm1] | |
36 | ||
37 | [graph tgamma] | |
38 | ||
39 | [optional_policy] | |
40 | ||
41 | There are effectively two versions of the [@http://en.wikipedia.org/wiki/Gamma_function tgamma] | |
42 | function internally: a fully | |
43 | generic version that is slow, but reasonably accurate, and a much more | |
44 | efficient approximation that is used where the number of digits in the significand | |
45 | of T correspond to a certain __lanczos. In practice any built in | |
46 | floating point type you will encounter has an appropriate __lanczos | |
47 | defined for it. It is also possible, given enough machine time, to generate | |
48 | further __lanczos's using the program libs/math/tools/lanczos_generator.cpp. | |
49 | ||
50 | The return type of this function is computed using the __arg_promotion_rules: | |
51 | the result is `double` when T is an integer type, and T otherwise. | |
52 | ||
53 | template <class T> | |
54 | ``__sf_result`` tgamma1pm1(T dz); | |
55 | ||
56 | template <class T, class ``__Policy``> | |
57 | ``__sf_result`` tgamma1pm1(T dz, const ``__Policy``&); | |
58 | ||
59 | Returns `tgamma(dz + 1) - 1`. Internally the implementation does not make | |
60 | use of the addition and subtraction implied by the definition, leading to | |
61 | accurate results even for very small `dz`. However, the implementation is | |
62 | capped to either 35 digit accuracy, or to the precision of the __lanczos | |
63 | associated with type T, whichever is more accurate. | |
64 | ||
65 | The return type of this function is computed using the __arg_promotion_rules: | |
66 | the result is `double` when T is an integer type, and T otherwise. | |
67 | ||
68 | [optional_policy] | |
69 | ||
70 | [h4 Accuracy] | |
71 | ||
72 | The following table shows the peak errors (in units of epsilon) | |
73 | found on various platforms with various floating point types, | |
74 | along with comparisons to other common libraries. | |
75 | Unless otherwise specified any floating point type that is narrower | |
76 | than the one shown will have __zero_error. | |
77 | ||
78 | [table_tgamma] | |
79 | ||
80 | [table_tgamma1pm1] | |
81 | ||
82 | [h4 Testing] | |
83 | ||
84 | The gamma is relatively easy to test: factorials and half-integer factorials | |
85 | can be calculated exactly by other means and compared with the gamma function. | |
86 | In addition, some accuracy tests in known tricky areas were computed at high precision | |
87 | using the generic version of this function. | |
88 | ||
89 | The function `tgamma1pm1` is tested against values calculated very naively | |
90 | using the formula `tgamma(1+dz)-1` with a lanczos approximation accurate | |
91 | to around 100 decimal digits. | |
92 | ||
93 | [h4 Implementation] | |
94 | ||
95 | The generic version of the `tgamma` function is implemented Sterling's approximation | |
96 | for lgamma for large z: | |
97 | ||
98 | [equation gamma6] | |
99 | ||
100 | Following exponentiation, downward recursion is then used for small values of z. | |
101 | ||
102 | For types of known precision the __lanczos is used, a traits class | |
103 | `boost::math::lanczos::lanczos_traits` maps type T to an appropriate | |
104 | approximation. | |
105 | ||
106 | For z in the range -20 < z < 1 then recursion is used to shift to z > 1 via: | |
107 | ||
108 | [equation gamm3] | |
109 | ||
110 | For very small z, this helps to preserve the identity: | |
111 | ||
112 | [equation gamm4] | |
113 | ||
114 | For z < -20 the reflection formula: | |
115 | ||
116 | [equation gamm5] | |
117 | ||
118 | is used. Particular care has to be taken to evaluate the [^ z * sin([pi][space] * z)] part: | |
119 | a special routine is used to reduce z prior to multiplying by [pi][space] to ensure that the | |
120 | result in is the range [0, [pi]/2]. Without this an excessive amount of error occurs | |
121 | in this region (which is hard enough already, as the rate of change near a negative pole | |
122 | is /exceptionally/ high). | |
123 | ||
124 | Finally if the argument is a small integer then table lookup of the factorial | |
125 | is used. | |
126 | ||
127 | The function `tgamma1pm1` is implemented using rational approximations [jm_rationals] in the | |
128 | region `-0.5 < dz < 2`. These are the same approximations (and internal routines) | |
129 | that are used for __lgamma, and so aren't detailed further here. The result of | |
130 | the approximation is `log(tgamma(dz+1))` which can fed into __expm1 to give | |
131 | the desired result. Outside the range `-0.5 < dz < 2` then the naive formula | |
132 | `tgamma1pm1(dz) = tgamma(dz+1)-1` can be used directly. | |
133 | ||
134 | [endsect][/section:tgamma The Gamma Function] | |
135 | [/ | |
136 | Copyright 2006 John Maddock and Paul A. Bristow. | |
137 | Distributed under the Boost Software License, Version 1.0. | |
138 | (See accompanying file LICENSE_1_0.txt or copy at | |
139 | http://www.boost.org/LICENSE_1_0.txt). | |
140 | ] | |
141 |