1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
[section:polygamma Polygamma]
[h4 Synopsis]
``
#include <boost/math/special_functions/polygamma.hpp>
``
namespace boost{ namespace math{
template <class T>
``__sf_result`` polygamma(int n, T z);
template <class T, class ``__Policy``>
``__sf_result`` polygamma(int n, T z, const ``__Policy``&);
}} // namespaces
[h4 Description]
Returns the polygamma function of /x/. Polygamma is defined as the n'th
derivative of the digamma function:
[equation polygamma1]
The following graphs illustrate the behaviour of the function for odd and even order:
[graph polygamma2]
[graph polygamma3]
[optional_policy]
The return type of this function is computed using the __arg_pomotion_rules:
the result is of type `double` when T is an integer type, and type T otherwise.
[h4 Accuracy]
The following table shows the peak errors (in units of epsilon)
found on various platforms with various floating point types.
Unless otherwise specified any floating point type that is narrower
than the one shown will have __zero_error.
[table
[[Significand Size] [Platform and Compiler] [Small-medium positive arguments] [Small-medium negative x] ]
[[53] [Win32 Visual C++ 12] [Peak=5.0 Mean=1] [Peak=1200 Mean=65]]
[[64] [Win64 Mingw GCC] [Peak=16 Mean=3] [Peak=33 Mean=3] ]
[[113] [Win64 Mingw GCC __float128] [Peak=6.5 Mean=1][Peak=30 Mean=4] ]
]
As shown above, error rates are generally very acceptable for moderately sized
arguments. Error rates should stay low for exact inputs, however, please note that the
function becomes exceptionally sensitive to small changes in input for large n and negative x,
indeed for cases where ['n!] would overflow, the function changes directly from -[infin] to
+[infin] somewhere between each negative integer - ['these cases are not handled correctly].
[*For these reasons results should be treated with extreme caution when /n/ is large and x negative].
[h4 Testing]
Testing is against Mathematica generated spot values to 35 digit precision.
[h4 Implementation]
For x < 0 the following reflection formula is used:
[equation polygamma2]
The n'th derivative of ['cot(x)] is tabulated for small /n/, and for larger n
has the general form:
[equation polygamma3]
The coefficients of the cosine terms can be calculated iteratively starting
from ['C[sub 1,0] = -1] and then using
[equation polygamma7]
to generate coefficients for n+1.
Note that every other coefficient is zero, and therefore what we have are
even or odd polynomials depending on whether n is even or odd.
Once x is positive then we have two methods available to us, for small x
we use the series expansion:
[equation polygamma4]
Note that the evaluation of zeta functions at integer values is essentially a table lookup
as __zeta is optimized for those cases.
For large x we use the asymptotic expansion:
[equation polygamma5]
For x in-between the two extremes we use the relation:
[equation polygamma6]
to make x large enough for the asymptotic expansion to be used.
There are also two special cases:
[equation polygamma8]
[equation polygamma9]
[endsect][/section:polygamma The Polygamma Function]
[/
Copyright 2014 John Maddock and Paul A. Bristow.
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or copy at
http://www.boost.org/LICENSE_1_0.txt).
]
|
