...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/digamma.hpp>
namespace boost{ namespace math{ template <class T> calculatedresulttype digamma(T z); template <class T, class Policy> calculatedresulttype digamma(T z, const Policy&); }} // namespaces
Returns the digamma or psi function of x. Digamma is defined as the logarithmic derivative of the gamma function:
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.4. Error rates for digamma
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Sun compiler version 0x5150 
Microsoft Visual C++ version 14.1 


Digamma Function: Large Values 
Max = 0ε (Mean = 0ε) 
Max = 1.39ε (Mean = 0.413ε) 
Max = 1.39ε (Mean = 0.413ε) 
Max = 0.98ε (Mean = 0.369ε) 
Digamma Function: Near the Positive Root 
Max = 0.891ε (Mean = 0.0995ε) 
Max = 1.37ε (Mean = 0.477ε) 
Max = 1.31ε (Mean = 0.471ε) 
Max = 0.997ε (Mean = 0.527ε) 
Digamma Function: Near Zero 
Max = 0ε (Mean = 0ε) 
Max = 0.984ε (Mean = 0.361ε) 
Max = 0.984ε (Mean = 0.361ε) 
Max = 0.953ε (Mean = 0.337ε) 
Digamma Function: Negative Values 
Max = 0ε (Mean = 0ε) 
Max = 180ε (Mean = 13ε) 
Max = 180ε (Mean = 13ε) 
Max = 214ε (Mean = 16.1ε) 
Digamma Function: Values near 0 
Max = 0ε (Mean = 0ε) 
Max = 1ε (Mean = 0.592ε) 
Max = 1ε (Mean = 0.592ε) 
Max = 0ε (Mean = 0ε) 
Digamma Function: Integer arguments 
Max = 0.992ε (Mean = 0.215ε) 
Max = 0.888ε (Mean = 0.403ε) 
Max = 0.888ε (Mean = 0.403ε) 
Max = 0.992ε (Mean = 0.452ε) 
Digamma Function: Half integer arguments 
Max = 0ε (Mean = 0ε) 
Max = 0.906ε (Mean = 0.409ε) 
Max = 0.906ε (Mean = 0.409ε) 
Max = 0.78ε (Mean = 0.314ε) 
As shown above, error rates for positive arguments are generally very low. For negative arguments there are an infinite number of irrational roots: relative errors very close to these can be arbitrarily large, although absolute error will remain very low.
The following error plot are based on an exhaustive search of the functions
domain, MSVC15.5 at double
precision, and GCC7.1/Ubuntu for long
double
and __float128
.
There are two sets of tests: spot values are computed using the online calculator at functions.wolfram.com, while random test values are generated using the highprecision reference implementation (a differentiated Lanczos approximation see below).
The implementation is divided up into the following domains:
For Negative arguments the reflection formula:
digamma(1x) = digamma(x) + pi/tan(pi*x);
is used to make x positive.
For arguments in the range [0,1] the recurrence relation:
digamma(x) = digamma(x+1)  1/x
is used to shift the evaluation to [1,2].
For arguments in the range [1,2] a rational approximation devised by JM is used (see below).
For arguments in the range [2,BIG] the recurrence relation:
digamma(x+1) = digamma(x) + 1/x;
is used to shift the evaluation to the range [1,2].
For arguments > BIG the asymptotic expansion:
can be used. However, this expansion is divergent after a few terms: exactly
how many terms depends on the size of x. Therefore the
value of BIG must be chosen so that the series can be
truncated at a term that is too small to have any effect on the result when
evaluated at BIG. Choosing BIG=10 for up to 80bit reals,
and BIG=20 for 128bit reals allows the series to truncated after a suitably
small number of terms and evaluated as a polynomial in 1/(x*x)
.
The arbitrary precision version of this function uses recurrence relations until x > BIG, and then evaluation via the asymptotic expansion above. As special cases integer and half integer arguments are handled via:
The rational approximation devised by JM in the range [1,2] is derived as follows.
First a high precision approximation to digamma was constructed using a 60term differentiated Lanczos approximation, the form used is:
Where P(x) and Q(x) are the polynomials from the rational form of the Lanczos
sum, and P'(x) and Q'(x) are their first derivatives. The Lanzos part of
this approximation has a theoretical precision of ~100 decimal digits. However,
cancellation in the above sum will reduce that to around 99(1/y)
digits
if y is the result. This approximation was used to calculate
the positive root of digamma, and was found to agree with the value used
by Cody to 25 digits (See Math. Comp. 27, 123127 (1973) by Cody, Strecok
and Thacher) and with the value used by Morris to 35 digits (See TOMS Algorithm
708).
Likewise a few spot tests agreed with values calculated using functions.wolfram.com to >40 digits. That's sufficiently precise to insure that the approximation below is accurate to double precision. Achieving 128bit long double precision requires that the location of the root is known to ~70 digits, and it's not clear whether the value calculated by this method meets that requirement: the difficulty lies in independently verifying the value obtained.
The rational approximation devised by JM was optimised for absolute error using the form:
digamma(x) = (x  X0)(Y + R(x  1));
Where X0 is the positive root of digamma, Y is a constant, and R(x  1) is the rational approximation. Note that since X0 is irrational, we need twice as many digits in X0 as in x in order to avoid cancellation error during the subtraction (this assumes that x is an exact value, if it's not then all bets are off). That means that even when x is the value of the root rounded to the nearest representable value, the result of digamma(x) will not be zero.