Boost C++ Libraries 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-type polygamma(int n, T z);

template <class T, class Policy>
calculated-result-type polygamma(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.

Table 8.6. Error rates for polygamma

GNU C++ version 7.1.0

GNU C++ version 7.1.0
long double

Sun compiler version 0x5150
Sun Solaris
long double

Microsoft Visual C++ version 14.1

Mathematica Data

Max = 0.824ε (Mean = 0.0574ε)

(GSL 2.1: Max = 62.9ε (Mean = 12.8ε))
(Rmath 3.2.3: Max = 108ε (Mean = 15.2ε))

Max = 7.38ε (Mean = 1.84ε)

Max = 34.3ε (Mean = 7.65ε)

Max = 9.32ε (Mean = 1.95ε)

Mathematica Data - large arguments

Max = 0.998ε (Mean = 0.0592ε)

(GSL 2.1: Max = 244ε (Mean = 32.8ε) And other failures.)
(Rmath 3.2.3: Max = 1.71e+56ε (Mean = 1.01e+55ε) And other failures.)

Max = 2.23ε (Mean = 0.323ε)

Max = 11.1ε (Mean = 0.848ε)

Max = 150ε (Mean = 13.9ε)

Mathematica Data - negative arguments

Max = 0.516ε (Mean = 0.022ε)

(GSL 2.1: Max = 36.6ε (Mean = 3.04ε) And other failures.)
(Rmath 3.2.3: Max = 0ε (Mean = 0ε) And other failures.)

Max = 269ε (Mean = 87.7ε)

Max = 269ε (Mean = 88.4ε)

Max = 497ε (Mean = 129ε)

Mathematica Data - large negative arguments

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 1.79ε (Mean = 0.197ε) And other failures.)
(Rmath 3.2.3: Max = 0ε (Mean = 0ε) And other failures.)

Max = 155ε (Mean = 96.4ε)

Max = 155ε (Mean = 96.4ε)

Max = 162ε (Mean = 101ε)

Mathematica Data - small arguments

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 15.2ε (Mean = 5.03ε))
(Rmath 3.2.3: Max = 106ε (Mean = 20ε))

Max = 3.33ε (Mean = 0.75ε)

Max = 3.33ε (Mean = 0.75ε)

Max = 3ε (Mean = 0.496ε)

Mathematica Data - Large orders and other bug cases

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 151ε (Mean = 39.3ε) And other failures.)
(Rmath 3.2.3: Max = +INFε (Mean = +INFε) And other failures.)

Max = 54.5ε (Mean = 13.3ε)

Max = 145ε (Mean = 55.9ε)

Max = 200ε (Mean = 57.2ε)

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 C1,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: