1[section:digamma Digamma] 2 3[h4 Synopsis] 4 5`` 6#include <boost/math/special_functions/digamma.hpp> 7`` 8 9 namespace boost{ namespace math{ 10 11 template <class T> 12 ``__sf_result`` digamma(T z); 13 14 template <class T, class ``__Policy``> 15 ``__sf_result`` digamma(T z, const ``__Policy``&); 16 17 }} // namespaces 18 19[h4 Description] 20 21Returns the digamma or psi function of /x/. Digamma is defined as the logarithmic 22derivative of the gamma function: 23 24[equation digamma1] 25 26[graph digamma] 27 28[optional_policy] 29 30The return type of this function is computed using the __arg_promotion_rules: 31the result is of type `double` when T is an integer type, and type T otherwise. 32 33[h4 Accuracy] 34 35The following table shows the peak errors (in units of epsilon) 36found on various platforms with various floating point types. 37Unless otherwise specified any floating point type that is narrower 38than the one shown will have __zero_error. 39 40[table_digamma] 41 42As shown above, error rates for positive arguments are generally very low. 43For negative arguments there are an infinite number of irrational roots: 44relative errors very close to these can be arbitrarily large, although 45absolute error will remain very low. 46 47The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision, 48and GCC-7.1/Ubuntu for `long double` and `__float128`. 49 50[graph digamma__double] 51 52[graph digamma__80_bit_long_double] 53 54[graph digamma____float128] 55 56[h4 Testing] 57 58There are two sets of tests: spot values are computed using 59the online calculator at functions.wolfram.com, while random test values 60are generated using the high-precision reference implementation (a 61differentiated __lanczos see below). 62 63[h4 Implementation] 64 65The implementation is divided up into the following domains: 66 67For Negative arguments the reflection formula: 68 69 digamma(1-x) = digamma(x) + pi/tan(pi*x); 70 71is used to make /x/ positive. 72 73For arguments in the range [0,1] the recurrence relation: 74 75 digamma(x) = digamma(x+1) - 1/x 76 77is used to shift the evaluation to [1,2]. 78 79For arguments in the range [1,2] a rational approximation [jm_rationals] is used (see below). 80 81For arguments in the range [2,BIG] the recurrence relation: 82 83 digamma(x+1) = digamma(x) + 1/x; 84 85is used to shift the evaluation to the range [1,2]. 86 87For arguments > BIG the asymptotic expansion: 88 89[equation digamma2] 90 91can be used. However, this expansion is divergent after a few terms: 92exactly how many terms depends on the size of /x/. Therefore the value 93of /BIG/ must be chosen so that the series can be truncated at a term 94that is too small to have any effect on the result when evaluated at /BIG/. 95Choosing BIG=10 for up to 80-bit reals, and BIG=20 for 128-bit reals allows 96the series to truncated after a suitably small number of terms and evaluated 97as a polynomial in `1/(x*x)`. 98 99The arbitrary precision version of this function uses recurrence relations until 100x > BIG, and then evaluation via the asymptotic expansion above. As special cases 101integer and half integer arguments are handled via: 102 103[equation digamma4] 104 105[equation digamma5] 106 107The rational approximation [jm_rationals] in the range [1,2] is derived as follows. 108 109First a high precision approximation to digamma was constructed using a 60-term 110differentiated __lanczos, the form used is: 111 112[equation digamma3] 113 114Where P(x) and Q(x) are the polynomials from the rational form of the Lanczos sum, 115and P'(x) and Q'(x) are their first derivatives. The Lanzos part of this 116approximation has a theoretical precision of ~100 decimal digits. However, 117cancellation in the above sum will reduce that to around `99-(1/y)` digits 118if /y/ is the result. This approximation was used to calculate the positive root 119of digamma, and was found to agree with the value used by 120Cody to 25 digits (See Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher) 121and with the value used by Morris to 35 digits (See TOMS Algorithm 708). 122 123Likewise a few spot tests agreed with values calculated using 124functions.wolfram.com to >40 digits. 125That's sufficiently precise to insure that the approximation below is 126accurate to double precision. Achieving 128-bit long double precision requires that 127the location of the root is known to ~70 digits, and it's not clear whether 128the value calculated by this method meets that requirement: the difficulty 129lies in independently verifying the value obtained. 130 131The rational approximation [jm_rationals] was optimised for absolute error using the form: 132 133 digamma(x) = (x - X0)(Y + R(x - 1)); 134 135Where X0 is the positive root of digamma, Y is a constant, and R(x - 1) is the 136rational approximation. Note that since X0 is irrational, we need twice as many 137digits in X0 as in x in order to avoid cancellation error during the subtraction 138(this assumes that /x/ is an exact value, if it's not then all bets are off). That 139means that even when x is the value of the root rounded to the nearest 140representable value, the result of digamma(x) ['[*will not be zero]]. 141 142 143[endsect] [/section:digamma The Digamma Function] 144 145[/ 146 Copyright 2006 John Maddock and Paul A. Bristow. 147 Distributed under the Boost Software License, Version 1.0. 148 (See accompanying file LICENSE_1_0.txt or copy at 149 http://www.boost.org/LICENSE_1_0.txt). 150] 151 152