libs/math/doc/sf/digamma.qbk

   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
  21 Returns the digamma or psi function of /x/. Digamma is defined as the logarithmic
  22 derivative of the gamma function:
  23
  24 [equation digamma1]
  25
  26 [graph digamma]
  27
  28 [optional_policy]
  29
  30 The return type of this function is computed using the __arg_promotion_rules:
  31 the result is of type `double` when T is an integer type, and type T otherwise.
  32
  33 [h4 Accuracy]
  34
  35 The following table shows the peak errors (in units of epsilon)
  36 found on various platforms with various floating point types.
  37 Unless otherwise specified any floating point type that is narrower
  38 than the one shown will have __zero_error.
  39
  40 [table_digamma]
  41
  42 As shown above, error rates for positive arguments are generally very low.
  43 For negative arguments there are an infinite number of irrational roots:
  44 relative errors very close to these can be arbitrarily large, although
  45 absolute error will remain very low.
  46
  47 The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision,
  48 and 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
  58 There are two sets of tests: spot values are computed using
  59 the online calculator at functions.wolfram.com, while random test values
  60 are generated using the high-precision reference implementation (a
  61 differentiated __lanczos see below).
  62
  63 [h4 Implementation]
  64
  65 The implementation is divided up into the following domains:
  66
  67 For Negative arguments the reflection formula:
  68
  69    digamma(1-x) = digamma(x) + pi/tan(pi*x);
  70
  71 is used to make /x/ positive.
  72
  73 For arguments in the range [0,1] the recurrence relation:
  74
  75    digamma(x) = digamma(x+1) - 1/x
  76
  77 is used to shift the evaluation to [1,2].
  78
  79 For arguments in the range [1,2] a rational approximation [jm_rationals] is used (see below).
  80
  81 For arguments in the range [2,BIG] the recurrence relation:
  82
  83    digamma(x+1) = digamma(x) + 1/x;
  84
  85 is used to shift the evaluation to the range [1,2].
  86
  87 For arguments > BIG the asymptotic expansion:
  88
  89 [equation digamma2]
  90
  91 can be used.  However, this expansion is divergent after a few terms:
  92 exactly how many terms depends on the size of /x/.  Therefore the value
  93 of /BIG/ must be chosen so that the series can be truncated at a term
  94 that is too small to have any effect on the result when evaluated at /BIG/.
  95 Choosing BIG=10 for up to 80-bit reals, and BIG=20 for 128-bit reals allows
  96 the series to truncated after a suitably small number of terms and evaluated
  97 as a polynomial in `1/(x*x)`.
  98
  99 The arbitrary precision version of this function uses recurrence relations until
 100 x > BIG, and then evaluation via the asymptotic expansion above.  As special cases
 101 integer and half integer arguments are handled via:
 102
 103 [equation digamma4]
 104
 105 [equation digamma5]
 106
 107 The rational approximation [jm_rationals] in the range [1,2] is derived as follows.
 108
 109 First a high precision approximation to digamma was constructed using a 60-term
 110 differentiated __lanczos, the form used is:
 111
 112 [equation digamma3]
 113
 114 Where P(x) and Q(x) are the polynomials from the rational form of the Lanczos sum,
 115 and P'(x) and Q'(x) are their first derivatives.  The Lanzos part of this
 116 approximation has a theoretical precision of ~100 decimal digits.  However,
 117 cancellation in the above sum will reduce that to around `99-(1/y)` digits
 118 if /y/ is the result.  This approximation was used to calculate the positive root
 119 of digamma, and was found to agree with the value used by
 120 Cody to 25 digits (See Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher)
 121 and with the value used by Morris to 35 digits (See TOMS Algorithm 708).
 122
 123 Likewise a few spot tests agreed with values calculated using
 124 functions.wolfram.com to >40 digits.
 125 That's sufficiently precise to insure that the approximation below is
 126 accurate to double precision.  Achieving 128-bit long double precision requires that
 127 the location of the root is known to ~70 digits, and it's not clear whether
 128 the value calculated by this method meets that requirement: the difficulty
 129 lies in independently verifying the value obtained.
 130
 131 The rational approximation [jm_rationals] was optimised for absolute error using the form:
 132
 133    digamma(x) = (x - X0)(Y + R(x - 1));
 134
 135 Where X0 is the positive root of digamma, Y is a constant, and R(x - 1) is the
 136 rational approximation.  Note that since X0 is irrational, we need twice as many
 137 digits 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
 139 means that even when x is the value of the root rounded to the nearest
 140 representable 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