Synopsis

#include <boost/math/special_functions/polygamma.hpp>

namespace boost{ namespace math{

template <class T>
calculated-result-type polygamma(int n, T z);

template <class T, class Policy>
calculated-result-type polygamma(int n, T z, const Policy&);

}} // namespaces

Description

Returns the polygamma function of x. Polygamma is defined as the n'th derivative of the digamma function:

The following graphs illustrate the behaviour of the function for odd and even order:

The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.

The return type of this function is computed using the result type calculation rules: the result is of type double when T is an integer type, and type T otherwise.

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 effectively zero error.

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 -∞ to +∞ 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.

Testing

Testing is against Mathematica generated spot values to 35 digit precision.

Implementation

For x < 0 the following reflection formula is used:

The n'th derivative of cot(x) is tabulated for small n, and for larger n has the general form:

The coefficients of the cosine terms can be calculated iteratively starting from C_1,0 = -1 and then using

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:

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:

For x in-between the two extremes we use the relation:

to make x large enough for the asymptotic expansion to be used.

There are also two special cases: