]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | [section:digamma Digamma] |
2 | ||
3 | [h4 Synopsis] | |
4 | ||
5 | `` | |
6 | #include <boost/math/special_functions/digamma.hpp> | |
7 | `` | |
8 | ||
9 | namespace boost{ namespace math{ | |
10 | ||
11 | template <class T> | |
12 | ``__sf_result`` digamma(T z); | |
13 | ||
14 | template <class T, class ``__Policy``> | |
15 | ``__sf_result`` digamma(T z, const ``__Policy``&); | |
16 | ||
17 | }} // namespaces | |
18 | ||
19 | [h4 Description] | |
20 | ||
21 | Returns the digamma or psi function of /x/. Digamma is defined as the logarithmic | |
22 | derivative of the gamma function: | |
23 | ||
24 | [equation digamma1] | |
25 | ||
26 | [graph digamma] | |
27 | ||
28 | [optional_policy] | |
29 | ||
30 | The return type of this function is computed using the __arg_promotion_rules: | |
31 | the result is of type `double` when T is an integer type, and type T otherwise. | |
32 | ||
33 | [h4 Accuracy] | |
34 | ||
35 | The following table shows the peak errors (in units of epsilon) | |
36 | found on various platforms with various floating point types. | |
37 | Unless otherwise specified any floating point type that is narrower | |
38 | than the one shown will have __zero_error. | |
39 | ||
40 | [table_digamma] | |
41 | ||
42 | As shown above, error rates for positive arguments are generally very low. | |
43 | For negative arguments there are an infinite number of irrational roots: | |
44 | relative errors very close to these can be arbitrarily large, although | |
45 | absolute error will remain very low. | |
46 | ||
47 | [h4 Testing] | |
48 | ||
49 | There are two sets of tests: spot values are computed using | |
50 | the online calculator at functions.wolfram.com, while random test values | |
51 | are generated using the high-precision reference implementation (a | |
52 | differentiated __lanczos see below). | |
53 | ||
54 | [h4 Implementation] | |
55 | ||
56 | The implementation is divided up into the following domains: | |
57 | ||
58 | For Negative arguments the reflection formula: | |
59 | ||
60 | digamma(1-x) = digamma(x) + pi/tan(pi*x); | |
61 | ||
62 | is used to make /x/ positive. | |
63 | ||
64 | For arguments in the range [0,1] the recurrence relation: | |
65 | ||
66 | digamma(x) = digamma(x+1) - 1/x | |
67 | ||
68 | is used to shift the evaluation to [1,2]. | |
69 | ||
70 | For arguments in the range [1,2] a rational approximation [jm_rationals] is used (see below). | |
71 | ||
72 | For arguments in the range [2,BIG] the recurrence relation: | |
73 | ||
74 | digamma(x+1) = digamma(x) + 1/x; | |
75 | ||
76 | is used to shift the evaluation to the range [1,2]. | |
77 | ||
78 | For arguments > BIG the asymptotic expansion: | |
79 | ||
80 | [equation digamma2] | |
81 | ||
82 | can be used. However, this expansion is divergent after a few terms: | |
83 | exactly how many terms depends on the size of /x/. Therefore the value | |
84 | of /BIG/ must be chosen so that the series can be truncated at a term | |
85 | that is too small to have any effect on the result when evaluated at /BIG/. | |
86 | Choosing BIG=10 for up to 80-bit reals, and BIG=20 for 128-bit reals allows | |
87 | the series to truncated after a suitably small number of terms and evaluated | |
88 | as a polynomial in `1/(x*x)`. | |
89 | ||
90 | The arbitrary precision version of this function uses recurrence relations until | |
91 | x > BIG, and then evaluation via the asymptotic expansion above. As special cases | |
92 | integer and half integer arguments are handled via: | |
93 | ||
94 | [equation digamma4] | |
95 | ||
96 | [equation digamma5] | |
97 | ||
98 | The rational approximation [jm_rationals] in the range [1,2] is derived as follows. | |
99 | ||
100 | First a high precision approximation to digamma was constructed using a 60-term | |
101 | differentiated __lanczos, the form used is: | |
102 | ||
103 | [equation digamma3] | |
104 | ||
105 | Where P(x) and Q(x) are the polynomials from the rational form of the Lanczos sum, | |
106 | and P'(x) and Q'(x) are their first derivatives. The Lanzos part of this | |
107 | approximation has a theoretical precision of ~100 decimal digits. However, | |
108 | cancellation in the above sum will reduce that to around `99-(1/y)` digits | |
109 | if /y/ is the result. This approximation was used to calculate the positive root | |
110 | of digamma, and was found to agree with the value used by | |
111 | Cody to 25 digits (See Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher) | |
112 | and with the value used by Morris to 35 digits (See TOMS Algorithm 708). | |
113 | ||
114 | Likewise a few spot tests agreed with values calculated using | |
115 | functions.wolfram.com to >40 digits. | |
116 | That's sufficiently precise to insure that the approximation below is | |
117 | accurate to double precision. Achieving 128-bit long double precision requires that | |
118 | the location of the root is known to ~70 digits, and it's not clear whether | |
119 | the value calculated by this method meets that requirement: the difficulty | |
120 | lies in independently verifying the value obtained. | |
121 | ||
122 | The rational approximation [jm_rationals] was optimised for absolute error using the form: | |
123 | ||
124 | digamma(x) = (x - X0)(Y + R(x - 1)); | |
125 | ||
126 | Where X0 is the positive root of digamma, Y is a constant, and R(x - 1) is the | |
127 | rational approximation. Note that since X0 is irrational, we need twice as many | |
128 | digits in X0 as in x in order to avoid cancellation error during the subtraction | |
129 | (this assumes that /x/ is an exact value, if it's not then all bets are off). That | |
130 | means that even when x is the value of the root rounded to the nearest | |
131 | representable value, the result of digamma(x) ['[*will not be zero]]. | |
132 | ||
133 | ||
134 | [endsect][/section:digamma The Digamma Function] | |
135 | ||
136 | [/ | |
137 | Copyright 2006 John Maddock and Paul A. Bristow. | |
138 | Distributed under the Boost Software License, Version 1.0. | |
139 | (See accompanying file LICENSE_1_0.txt or copy at | |
140 | http://www.boost.org/LICENSE_1_0.txt). | |
141 | ] | |
142 |