...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> calculatedresulttype polygamma(int n, T z); template <class T, class Policy> calculatedresulttype 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 6.6. Error rates for polygamma
Microsoft Visual C++ version 12.0 
GNU C++ version 5.1.0 
GNU C++ version 5.1.0 
Sun compiler version 0x5130 


Mathematica Data 
Max = 6.34ε (Mean = 1.53ε) 
Max = 0.824ε (Mean = 0.0574ε) 
Max = 7.38ε (Mean = 1.84ε) 
Max = 18.3ε (Mean = 4.16ε) 
Mathematica Data  large arguments 
Max = 150ε (Mean = 15.1ε) 
Max = 0.998ε (Mean = 0.0592ε) 
Max = 2.23ε (Mean = 0.323ε) 
Max = 2.35ε (Mean = 0.34ε) 
Mathematica Data  negative arguments 
Max = 497ε (Mean = 129ε) 
Max = 0.516ε (Mean = 0.022ε) 
Max = 269ε (Mean = 87.7ε) 
Max = 269ε (Mean = 87.9ε) 
Mathematica Data  large negative arguments 
Max = 162ε (Mean = 101ε) 
Max = 0ε (Mean = 0ε) 
Max = 155ε (Mean = 96.4ε) 
Max = 155ε (Mean = 96.4ε) 
Mathematica Data  small arguments 
Max = 3ε (Mean = 0.496ε) 
Max = 0ε (Mean = 0ε) 
Max = 3.33ε (Mean = 0.75ε) 
Max = 3.33ε (Mean = 0.75ε) 
Mathematica Data  Large orders and other bug cases 
Max = 200ε (Mean = 57.2ε) 
Max = 0ε (Mean = 0ε) 
Max = 54.5ε (Mean = 13.3ε) 
Max = 90.1ε (Mean = 30.6ε) 
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 inbetween 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: