1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Digamma</title> 5<link rel="stylesheet" href="../../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../sf_gamma.html" title="Gamma Functions"> 9<link rel="prev" href="lgamma.html" title="Log Gamma"> 10<link rel="next" href="trigamma.html" title="Trigamma"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="lgamma.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../sf_gamma.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="trigamma.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="math_toolkit.sf_gamma.digamma"></a><a class="link" href="digamma.html" title="Digamma">Digamma</a> 28</h3></div></div></div> 29<h5> 30<a name="math_toolkit.sf_gamma.digamma.h0"></a> 31 <span class="phrase"><a name="math_toolkit.sf_gamma.digamma.synopsis"></a></span><a class="link" href="digamma.html#math_toolkit.sf_gamma.digamma.synopsis">Synopsis</a> 32 </h5> 33<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">digamma</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 34</pre> 35<pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span><span class="special">{</span> 36 37<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 38<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">digamma</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">);</span> 39 40<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 41<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">digamma</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> 42 43<span class="special">}}</span> <span class="comment">// namespaces</span> 44</pre> 45<h5> 46<a name="math_toolkit.sf_gamma.digamma.h1"></a> 47 <span class="phrase"><a name="math_toolkit.sf_gamma.digamma.description"></a></span><a class="link" href="digamma.html#math_toolkit.sf_gamma.digamma.description">Description</a> 48 </h5> 49<p> 50 Returns the digamma or psi function of <span class="emphasis"><em>x</em></span>. Digamma is 51 defined as the logarithmic derivative of the gamma function: 52 </p> 53<div class="blockquote"><blockquote class="blockquote"><p> 54 <span class="inlinemediaobject"><img src="../../../equations/digamma1.svg"></span> 55 56 </p></blockquote></div> 57<div class="blockquote"><blockquote class="blockquote"><p> 58 <span class="inlinemediaobject"><img src="../../../graphs/digamma.svg" align="middle"></span> 59 60 </p></blockquote></div> 61<p> 62 The final <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a> argument is optional and can 63 be used to control the behaviour of the function: how it handles errors, 64 what level of precision to use etc. Refer to the <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">policy 65 documentation for more details</a>. 66 </p> 67<p> 68 The return type of this function is computed using the <a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>result 69 type calculation rules</em></span></a>: the result is of type <code class="computeroutput"><span class="keyword">double</span></code> when T is an integer type, and type 70 T otherwise. 71 </p> 72<h5> 73<a name="math_toolkit.sf_gamma.digamma.h2"></a> 74 <span class="phrase"><a name="math_toolkit.sf_gamma.digamma.accuracy"></a></span><a class="link" href="digamma.html#math_toolkit.sf_gamma.digamma.accuracy">Accuracy</a> 75 </h5> 76<p> 77 The following table shows the peak errors (in units of epsilon) found on 78 various platforms with various floating point types. Unless otherwise specified 79 any floating point type that is narrower than the one shown will have <a class="link" href="../relative_error.html#math_toolkit.relative_error.zero_error">effectively zero error</a>. 80 </p> 81<div class="table"> 82<a name="math_toolkit.sf_gamma.digamma.table_digamma"></a><p class="title"><b>Table 8.4. Error rates for digamma</b></p> 83<div class="table-contents"><table class="table" summary="Error rates for digamma"> 84<colgroup> 85<col> 86<col> 87<col> 88<col> 89<col> 90</colgroup> 91<thead><tr> 92<th> 93 </th> 94<th> 95 <p> 96 GNU C++ version 7.1.0<br> linux<br> double 97 </p> 98 </th> 99<th> 100 <p> 101 GNU C++ version 7.1.0<br> linux<br> long double 102 </p> 103 </th> 104<th> 105 <p> 106 Sun compiler version 0x5150<br> Sun Solaris<br> long double 107 </p> 108 </th> 109<th> 110 <p> 111 Microsoft Visual C++ version 14.1<br> Win32<br> double 112 </p> 113 </th> 114</tr></thead> 115<tbody> 116<tr> 117<td> 118 <p> 119 Digamma Function: Large Values 120 </p> 121 </td> 122<td> 123 <p> 124 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 125 2.1:</em></span> Max = 1.84ε (Mean = 0.71ε))<br> (<span class="emphasis"><em>Rmath 126 3.2.3:</em></span> Max = 1.18ε (Mean = 0.331ε)) 127 </p> 128 </td> 129<td> 130 <p> 131 <span class="blue">Max = 1.39ε (Mean = 0.413ε)</span> 132 </p> 133 </td> 134<td> 135 <p> 136 <span class="blue">Max = 1.39ε (Mean = 0.413ε)</span> 137 </p> 138 </td> 139<td> 140 <p> 141 <span class="blue">Max = 0.98ε (Mean = 0.369ε)</span> 142 </p> 143 </td> 144</tr> 145<tr> 146<td> 147 <p> 148 Digamma Function: Near the Positive Root 149 </p> 150 </td> 151<td> 152 <p> 153 <span class="blue">Max = 0.891ε (Mean = 0.0995ε)</span><br> 154 <br> (<span class="emphasis"><em>GSL 2.1:</em></span> Max = 135ε (Mean = 11.9ε))<br> 155 (<span class="emphasis"><em>Rmath 3.2.3:</em></span> Max = 2.02e+03ε (Mean = 256ε)) 156 </p> 157 </td> 158<td> 159 <p> 160 <span class="blue">Max = 1.37ε (Mean = 0.477ε)</span> 161 </p> 162 </td> 163<td> 164 <p> 165 <span class="blue">Max = 1.31ε (Mean = 0.471ε)</span> 166 </p> 167 </td> 168<td> 169 <p> 170 <span class="blue">Max = 0.997ε (Mean = 0.527ε)</span> 171 </p> 172 </td> 173</tr> 174<tr> 175<td> 176 <p> 177 Digamma Function: Near Zero 178 </p> 179 </td> 180<td> 181 <p> 182 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 183 2.1:</em></span> Max = 0.953ε (Mean = 0.348ε))<br> (<span class="emphasis"><em>Rmath 184 3.2.3:</em></span> Max = 1.17ε (Mean = 0.564ε)) 185 </p> 186 </td> 187<td> 188 <p> 189 <span class="blue">Max = 0.984ε (Mean = 0.361ε)</span> 190 </p> 191 </td> 192<td> 193 <p> 194 <span class="blue">Max = 0.984ε (Mean = 0.361ε)</span> 195 </p> 196 </td> 197<td> 198 <p> 199 <span class="blue">Max = 0.953ε (Mean = 0.337ε)</span> 200 </p> 201 </td> 202</tr> 203<tr> 204<td> 205 <p> 206 Digamma Function: Negative Values 207 </p> 208 </td> 209<td> 210 <p> 211 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 212 2.1:</em></span> Max = 4.56e+04ε (Mean = 3.91e+03ε))<br> (<span class="emphasis"><em>Rmath 213 3.2.3:</em></span> Max = 4.6e+04ε (Mean = 3.94e+03ε)) 214 </p> 215 </td> 216<td> 217 <p> 218 <span class="blue">Max = 180ε (Mean = 13ε)</span> 219 </p> 220 </td> 221<td> 222 <p> 223 <span class="blue">Max = 180ε (Mean = 13ε)</span> 224 </p> 225 </td> 226<td> 227 <p> 228 <span class="blue">Max = 214ε (Mean = 16.1ε)</span> 229 </p> 230 </td> 231</tr> 232<tr> 233<td> 234 <p> 235 Digamma Function: Values near 0 236 </p> 237 </td> 238<td> 239 <p> 240 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 241 2.1:</em></span> Max = 0.866ε (Mean = 0.387ε))<br> (<span class="emphasis"><em>Rmath 242 3.2.3:</em></span> Max = 3.58e+05ε (Mean = 1.6e+05ε)) 243 </p> 244 </td> 245<td> 246 <p> 247 <span class="blue">Max = 1ε (Mean = 0.592ε)</span> 248 </p> 249 </td> 250<td> 251 <p> 252 <span class="blue">Max = 1ε (Mean = 0.592ε)</span> 253 </p> 254 </td> 255<td> 256 <p> 257 <span class="blue">Max = 0ε (Mean = 0ε)</span> 258 </p> 259 </td> 260</tr> 261<tr> 262<td> 263 <p> 264 Digamma Function: Integer arguments 265 </p> 266 </td> 267<td> 268 <p> 269 <span class="blue">Max = 0.992ε (Mean = 0.215ε)</span><br> <br> 270 (<span class="emphasis"><em>GSL 2.1:</em></span> Max = 1.18ε (Mean = 0.607ε))<br> 271 (<span class="emphasis"><em>Rmath 3.2.3:</em></span> Max = 4.33ε (Mean = 0.982ε)) 272 </p> 273 </td> 274<td> 275 <p> 276 <span class="blue">Max = 0.888ε (Mean = 0.403ε)</span> 277 </p> 278 </td> 279<td> 280 <p> 281 <span class="blue">Max = 0.888ε (Mean = 0.403ε)</span> 282 </p> 283 </td> 284<td> 285 <p> 286 <span class="blue">Max = 0.992ε (Mean = 0.452ε)</span> 287 </p> 288 </td> 289</tr> 290<tr> 291<td> 292 <p> 293 Digamma Function: Half integer arguments 294 </p> 295 </td> 296<td> 297 <p> 298 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 299 2.1:</em></span> Max = 1.09ε (Mean = 0.531ε))<br> (<span class="emphasis"><em>Rmath 300 3.2.3:</em></span> Max = 46.2ε (Mean = 7.24ε)) 301 </p> 302 </td> 303<td> 304 <p> 305 <span class="blue">Max = 0.906ε (Mean = 0.409ε)</span> 306 </p> 307 </td> 308<td> 309 <p> 310 <span class="blue">Max = 0.906ε (Mean = 0.409ε)</span> 311 </p> 312 </td> 313<td> 314 <p> 315 <span class="blue">Max = 0.78ε (Mean = 0.314ε)</span> 316 </p> 317 </td> 318</tr> 319</tbody> 320</table></div> 321</div> 322<br class="table-break"><p> 323 As shown above, error rates for positive arguments are generally very low. 324 For negative arguments there are an infinite number of irrational roots: 325 relative errors very close to these can be arbitrarily large, although absolute 326 error will remain very low. 327 </p> 328<p> 329 The following error plot are based on an exhaustive search of the functions 330 domain, MSVC-15.5 at <code class="computeroutput"><span class="keyword">double</span></code> 331 precision, and GCC-7.1/Ubuntu for <code class="computeroutput"><span class="keyword">long</span> 332 <span class="keyword">double</span></code> and <code class="computeroutput"><span class="identifier">__float128</span></code>. 333 </p> 334<div class="blockquote"><blockquote class="blockquote"><p> 335 <span class="inlinemediaobject"><img src="../../../graphs/digamma__double.svg" align="middle"></span> 336 337 </p></blockquote></div> 338<div class="blockquote"><blockquote class="blockquote"><p> 339 <span class="inlinemediaobject"><img src="../../../graphs/digamma__80_bit_long_double.svg" align="middle"></span> 340 341 </p></blockquote></div> 342<div class="blockquote"><blockquote class="blockquote"><p> 343 <span class="inlinemediaobject"><img src="../../../graphs/digamma____float128.svg" align="middle"></span> 344 345 </p></blockquote></div> 346<h5> 347<a name="math_toolkit.sf_gamma.digamma.h3"></a> 348 <span class="phrase"><a name="math_toolkit.sf_gamma.digamma.testing"></a></span><a class="link" href="digamma.html#math_toolkit.sf_gamma.digamma.testing">Testing</a> 349 </h5> 350<p> 351 There are two sets of tests: spot values are computed using the online calculator 352 at functions.wolfram.com, while random test values are generated using the 353 high-precision reference implementation (a differentiated <a class="link" href="../lanczos.html" title="The Lanczos Approximation">Lanczos 354 approximation</a> see below). 355 </p> 356<h5> 357<a name="math_toolkit.sf_gamma.digamma.h4"></a> 358 <span class="phrase"><a name="math_toolkit.sf_gamma.digamma.implementation"></a></span><a class="link" href="digamma.html#math_toolkit.sf_gamma.digamma.implementation">Implementation</a> 359 </h5> 360<p> 361 The implementation is divided up into the following domains: 362 </p> 363<p> 364 For Negative arguments the reflection formula: 365 </p> 366<pre class="programlisting"><span class="identifier">digamma</span><span class="special">(</span><span class="number">1</span><span class="special">-</span><span class="identifier">x</span><span class="special">)</span> <span class="special">=</span> <span class="identifier">digamma</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special">+</span> <span class="identifier">pi</span><span class="special">/</span><span class="identifier">tan</span><span class="special">(</span><span class="identifier">pi</span><span class="special">*</span><span class="identifier">x</span><span class="special">);</span> 367</pre> 368<p> 369 is used to make <span class="emphasis"><em>x</em></span> positive. 370 </p> 371<p> 372 For arguments in the range [0,1] the recurrence relation: 373 </p> 374<pre class="programlisting"><span class="identifier">digamma</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special">=</span> <span class="identifier">digamma</span><span class="special">(</span><span class="identifier">x</span><span class="special">+</span><span class="number">1</span><span class="special">)</span> <span class="special">-</span> <span class="number">1</span><span class="special">/</span><span class="identifier">x</span> 375</pre> 376<p> 377 is used to shift the evaluation to [1,2]. 378 </p> 379<p> 380 For arguments in the range [1,2] a rational approximation <a class="link" href="../sf_implementation.html#math_toolkit.sf_implementation.rational_approximations_used">devised 381 by JM</a> is used (see below). 382 </p> 383<p> 384 For arguments in the range [2,BIG] the recurrence relation: 385 </p> 386<pre class="programlisting"><span class="identifier">digamma</span><span class="special">(</span><span class="identifier">x</span><span class="special">+</span><span class="number">1</span><span class="special">)</span> <span class="special">=</span> <span class="identifier">digamma</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special">+</span> <span class="number">1</span><span class="special">/</span><span class="identifier">x</span><span class="special">;</span> 387</pre> 388<p> 389 is used to shift the evaluation to the range [1,2]. 390 </p> 391<p> 392 For arguments > BIG the asymptotic expansion: 393 </p> 394<div class="blockquote"><blockquote class="blockquote"><p> 395 <span class="inlinemediaobject"><img src="../../../equations/digamma2.svg"></span> 396 397 </p></blockquote></div> 398<p> 399 can be used. However, this expansion is divergent after a few terms: exactly 400 how many terms depends on the size of <span class="emphasis"><em>x</em></span>. Therefore the 401 value of <span class="emphasis"><em>BIG</em></span> must be chosen so that the series can be 402 truncated at a term that is too small to have any effect on the result when 403 evaluated at <span class="emphasis"><em>BIG</em></span>. Choosing BIG=10 for up to 80-bit reals, 404 and BIG=20 for 128-bit reals allows the series to truncated after a suitably 405 small number of terms and evaluated as a polynomial in <code class="computeroutput"><span class="number">1</span><span class="special">/(</span><span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span><span class="special">)</span></code>. 406 </p> 407<p> 408 The arbitrary precision version of this function uses recurrence relations 409 until x > BIG, and then evaluation via the asymptotic expansion above. 410 As special cases integer and half integer arguments are handled via: 411 </p> 412<div class="blockquote"><blockquote class="blockquote"><p> 413 <span class="inlinemediaobject"><img src="../../../equations/digamma4.svg"></span> 414 415 </p></blockquote></div> 416<div class="blockquote"><blockquote class="blockquote"><p> 417 <span class="inlinemediaobject"><img src="../../../equations/digamma5.svg"></span> 418 419 </p></blockquote></div> 420<p> 421 The rational approximation <a class="link" href="../sf_implementation.html#math_toolkit.sf_implementation.rational_approximations_used">devised 422 by JM</a> in the range [1,2] is derived as follows. 423 </p> 424<p> 425 First a high precision approximation to digamma was constructed using a 60-term 426 differentiated <a class="link" href="../lanczos.html" title="The Lanczos Approximation">Lanczos approximation</a>, 427 the form used is: 428 </p> 429<div class="blockquote"><blockquote class="blockquote"><p> 430 <span class="inlinemediaobject"><img src="../../../equations/digamma3.svg"></span> 431 432 </p></blockquote></div> 433<p> 434 Where P(x) and Q(x) are the polynomials from the rational form of the Lanczos 435 sum, and P'(x) and Q'(x) are their first derivatives. The Lanzos part of 436 this approximation has a theoretical precision of ~100 decimal digits. However, 437 cancellation in the above sum will reduce that to around <code class="computeroutput"><span class="number">99</span><span class="special">-(</span><span class="number">1</span><span class="special">/</span><span class="identifier">y</span><span class="special">)</span></code> digits 438 if <span class="emphasis"><em>y</em></span> is the result. This approximation was used to calculate 439 the positive root of digamma, and was found to agree with the value used 440 by Cody to 25 digits (See Math. Comp. 27, 123-127 (1973) by Cody, Strecok 441 and Thacher) and with the value used by Morris to 35 digits (See TOMS Algorithm 442 708). 443 </p> 444<p> 445 Likewise a few spot tests agreed with values calculated using functions.wolfram.com 446 to >40 digits. That's sufficiently precise to insure that the approximation 447 below is accurate to double precision. Achieving 128-bit long double precision 448 requires that the location of the root is known to ~70 digits, and it's not 449 clear whether the value calculated by this method meets that requirement: 450 the difficulty lies in independently verifying the value obtained. 451 </p> 452<p> 453 The rational approximation <a class="link" href="../sf_implementation.html#math_toolkit.sf_implementation.rational_approximations_used">devised 454 by JM</a> was optimised for absolute error using the form: 455 </p> 456<pre class="programlisting"><span class="identifier">digamma</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">X0</span><span class="special">)(</span><span class="identifier">Y</span> <span class="special">+</span> <span class="identifier">R</span><span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="number">1</span><span class="special">));</span> 457</pre> 458<p> 459 Where X0 is the positive root of digamma, Y is a constant, and R(x - 1) is 460 the rational approximation. Note that since X0 is irrational, we need twice 461 as many digits in X0 as in x in order to avoid cancellation error during 462 the subtraction (this assumes that <span class="emphasis"><em>x</em></span> is an exact value, 463 if it's not then all bets are off). That means that even when x is the value 464 of the root rounded to the nearest representable value, the result of digamma(x) 465 <span class="emphasis"><em><span class="bold"><strong>will not be zero</strong></span></em></span>. 466 </p> 467</div> 468<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 469<td align="left"></td> 470<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 471 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 472 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 473 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 474 Daryle Walker and Xiaogang Zhang<p> 475 Distributed under the Boost Software License, Version 1.0. (See accompanying 476 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 477 </p> 478</div></td> 479</tr></table> 480<hr> 481<div class="spirit-nav"> 482<a accesskey="p" href="lgamma.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../sf_gamma.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="trigamma.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 483</div> 484</body> 485</html> 486
