...one of the most highly
regarded and expertly designed C++ library projects in the
world.

— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards

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

namespace boost{ namespace math{ template <class T>calculated-result-typepolygamma(int n, T z); template <class T, class Policy>calculated-result-typepolygamma(int n, T z, const Policy&); }} // namespaces

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.

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.

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 -∞ 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 is against Mathematica generated spot values to 35 digit precision.

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: