<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Digamma</title>
<link rel="stylesheet" href="../../../../../../../../doc/src/boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.74.0">
<link rel="home" href="../../../index.html" title="Math Toolkit">
<link rel="up" href="../sf_gamma.html" title="Gamma Functions">
<link rel="prev" href="lgamma.html" title="Log Gamma">
<link rel="next" href="gamma_ratios.html" title="Ratios of Gamma Functions">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<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="gamma_ratios.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section" lang="en">
<div class="titlepage"><div><div><h4 class="title">
<a name="math_toolkit.special.sf_gamma.digamma"></a><a class="link" href="digamma.html" title="Digamma"> Digamma</a>
</h4></div></div></div>
<a name="math_toolkit.special.sf_gamma.digamma.synopsis"></a><h5>
<a name="id1076683"></a>
<a class="link" href="digamma.html#math_toolkit.special.sf_gamma.digamma.synopsis">Synopsis</a>
</h5>
<p>
</p>
<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>
</pre>
<p>
</p>
<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>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span>
<a class="link" href="../../main_overview/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">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="Policies">Policy</a><span class="special">></span>
<a class="link" href="../../main_overview/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="Policies">Policy</a><span class="special">&);</span>
<span class="special">}}</span> <span class="comment">// namespaces
</span></pre>
<a name="math_toolkit.special.sf_gamma.digamma.description"></a><h5>
<a name="id1076910"></a>
<a class="link" href="digamma.html#math_toolkit.special.sf_gamma.digamma.description">Description</a>
</h5>
<p>
Returns the digamma or psi function of <span class="emphasis"><em>x</em></span>. Digamma
is defined as the logarithmic derivative of the gamma function:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/digamma1.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../graphs/digamma.png" align="middle"></span>
</p>
<p>
</p>
<p>
The final <a class="link" href="../../policy.html" title="Policies">Policy</a> 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
<a class="link" href="../../policy.html" title="Policies">policy documentation for more details</a>.
</p>
<p>
</p>
<p>
There is no fully generic version of this function: all the implementations
are tuned to specific accuracy levels, the most precise of which delivers
34-digits of precision.
</p>
<p>
The return type of this function is computed using the <a class="link" href="../../main_overview/result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>result
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
T otherwise.
</p>
<a name="math_toolkit.special.sf_gamma.digamma.accuracy"></a><h5>
<a name="id1077018"></a>
<a class="link" href="digamma.html#math_toolkit.special.sf_gamma.digamma.accuracy">Accuracy</a>
</h5>
<p>
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
<a class="link" href="../../backgrounders/relative_error.html#zero_error">effectively zero error</a>.
</p>
<div class="informaltable"><table class="table">
<colgroup>
<col>
<col>
<col>
<col>
<col>
<col>
</colgroup>
<thead><tr>
<th>
<p>
Significand Size
</p>
</th>
<th>
<p>
Platform and Compiler
</p>
</th>
<th>
<p>
Random Positive Values
</p>
</th>
<th>
<p>
Values Near The Positive Root
</p>
</th>
<th>
<p>
Values Near Zero
</p>
</th>
<th>
<p>
Negative Values
</p>
</th>
</tr></thead>
<tbody>
<tr>
<td>
<p>
53
</p>
</td>
<td>
<p>
Win32 Visual C++ 8
</p>
</td>
<td>
<p>
Peak=0.98 Mean=0.36
</p>
</td>
<td>
<p>
Peak=0.99 Mean=0.5
</p>
</td>
<td>
<p>
Peak=0.95 Mean=0.5
</p>
</td>
<td>
<p>
Peak=214 Mean=16
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
Linux IA32 / GCC
</p>
</td>
<td>
<p>
Peak=1.4 Mean=0.4
</p>
</td>
<td>
<p>
Peak=1.3 Mean=0.45
</p>
</td>
<td>
<p>
Peak=0.98 Mean=0.35
</p>
</td>
<td>
<p>
Peak=180 Mean=13
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
Linux IA64 / GCC
</p>
</td>
<td>
<p>
Peak=0.92 Mean=0.4
</p>
</td>
<td>
<p>
Peak=1.3 Mean=0.45
</p>
</td>
<td>
<p>
Peak=0.98 Mean=0.4
</p>
</td>
<td>
<p>
Peak=180 Mean=13
</p>
</td>
</tr>
<tr>
<td>
<p>
113
</p>
</td>
<td>
<p>
HPUX IA64, aCC A.06.06
</p>
</td>
<td>
<p>
Peak=0.9 Mean=0.4
</p>
</td>
<td>
<p>
Peak=1.1 Mean=0.5
</p>
</td>
<td>
<p>
Peak=0.99 Mean=0.4
</p>
</td>
<td>
<p>
Peak=64 Mean=6
</p>
</td>
</tr>
</tbody>
</table></div>
<p>
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.
</p>
<a name="math_toolkit.special.sf_gamma.digamma.testing"></a><h5>
<a name="id1077322"></a>
<a class="link" href="digamma.html#math_toolkit.special.sf_gamma.digamma.testing">Testing</a>
</h5>
<p>
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 high-precision reference implementation (a differentiated <a class="link" href="../../backgrounders/lanczos.html" title="The Lanczos Approximation">Lanczos approximation</a>
see below).
</p>
<a name="math_toolkit.special.sf_gamma.digamma.implementation"></a><h5>
<a name="id1077343"></a>
<a class="link" href="digamma.html#math_toolkit.special.sf_gamma.digamma.implementation">Implementation</a>
</h5>
<p>
The implementation is divided up into the following domains:
</p>
<p>
For Negative arguments the reflection formula:
</p>
<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>
</pre>
<p>
is used to make <span class="emphasis"><em>x</em></span> positive.
</p>
<p>
For arguments in the range [0,1] the recurrence relation:
</p>
<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>
</pre>
<p>
is used to shift the evaluation to [1,2].
</p>
<p>
For arguments in the range [1,2] a rational approximation <a class="link" href="../../backgrounders/implementation.html#math_toolkit.backgrounders.implementation.rational_approximations_used">devised
by JM</a> is used (see below).
</p>
<p>
For arguments in the range [2,BIG] the recurrence relation:
</p>
<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>
</pre>
<p>
is used to shift the evaluation to the range [1,2].
</p>
<p>
For arguments > BIG the asymptotic expansion:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/digamma2.png"></span>
</p>
<p>
can be used. However, this expansion is divergent after a few terms: exactly
how many terms depends on the size of <span class="emphasis"><em>x</em></span>. Therefore
the value of <span class="emphasis"><em>BIG</em></span> 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 <span class="emphasis"><em>BIG</em></span>. Choosing BIG=10 for
up to 80-bit reals, and BIG=20 for 128-bit reals allows the series to truncated
after a suitably 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>.
</p>
<p>
The rational approximation <a class="link" href="../../backgrounders/implementation.html#math_toolkit.backgrounders.implementation.rational_approximations_used">devised
by JM</a> in the range [1,2] is derived as follows.
</p>
<p>
First a high precision approximation to digamma was constructed using a
60-term differentiated <a class="link" href="../../backgrounders/lanczos.html" title="The Lanczos Approximation">Lanczos
approximation</a>, the form used is:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/digamma3.png"></span>
</p>
<p>
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 <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 if <span class="emphasis"><em>y</em></span> 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, 123-127 (1973)
by Cody, Strecok and Thacher) and with the value used by Morris to 35 digits
(See TOMS Algorithm 708).
</p>
<p>
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 128-bit 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.
</p>
<p>
The rational approximation <a class="link" href="../../backgrounders/implementation.html#math_toolkit.backgrounders.implementation.rational_approximations_used">devised
by JM</a> was optimised for absolute error using the form:
</p>
<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>
</pre>
<p>
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 <span class="emphasis"><em>x</em></span> 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) <span class="emphasis"><em><span class="bold"><strong>will not be zero</strong></span></em></span>.
</p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006 , 2007, 2008, 2009 John Maddock, Paul A. Bristow,
Hubert Holin, Xiaogang Zhang, Bruno Lalande, Johan Råde, Gautam Sewani
and Thijs van den Berg<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
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>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<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="gamma_ratios.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
