libs/math/doc/html/math_toolkit/sf_gamma/lgamma.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   4 <title>Log Gamma</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.11.0">
   8 <link rel="up" href="../sf_gamma.html" title="Gamma Functions">
   9 <link rel="prev" href="tgamma.html" title="Gamma">
  10 <link rel="next" href="digamma.html" title="Digamma">
  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="tgamma.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="digamma.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.lgamma"></a><a class="link" href="lgamma.html" title="Log Gamma">Log Gamma</a>
  28 </h3></div></div></div>
  29 <h5>
  30 <a name="math_toolkit.sf_gamma.lgamma.h0"></a>
  31         <span class="phrase"><a name="math_toolkit.sf_gamma.lgamma.synopsis"></a></span><a class="link" href="lgamma.html#math_toolkit.sf_gamma.lgamma.synopsis">Synopsis</a>
  32       </h5>
  33 <pre class="programlisting"><span class="preprocessor">#include</span> <span class="special">&lt;</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">gamma</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">&gt;</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">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</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">lgamma</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">&lt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&gt;</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">lgamma</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&amp;);</span>
  42
  43 <span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</span>
  44 <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">lgamma</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">int</span><span class="special">*</span> <span class="identifier">sign</span><span class="special">);</span>
  45
  46 <span class="keyword">template</span> <span class="special">&lt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&gt;</span>
  47 <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">lgamma</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">int</span><span class="special">*</span> <span class="identifier">sign</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../../policy.html" title="Chapter&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&amp;);</span>
  48
  49 <span class="special">}}</span> <span class="comment">// namespaces</span>
  50 </pre>
  51 <h5>
  52 <a name="math_toolkit.sf_gamma.lgamma.h1"></a>
  53         <span class="phrase"><a name="math_toolkit.sf_gamma.lgamma.description"></a></span><a class="link" href="lgamma.html#math_toolkit.sf_gamma.lgamma.description">Description</a>
  54       </h5>
  55 <p>
  56         The <a href="http://en.wikipedia.org/wiki/Gamma_function" target="_top">lgamma function</a>
  57         is defined by:
  58       </p>
  59 <div class="blockquote"><blockquote class="blockquote"><p>
  60           <span class="inlinemediaobject"><img src="../../../equations/lgamm1.svg"></span>
  61
  62         </p></blockquote></div>
  63 <p>
  64         The second form of the function takes a pointer to an integer, which if non-null
  65         is set on output to the sign of tgamma(z).
  66       </p>
  67 <p>
  68         The final <a class="link" href="../../policy.html" title="Chapter&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a> argument is optional and can
  69         be used to control the behaviour of the function: how it handles errors,
  70         what level of precision to use etc. Refer to the <a class="link" href="../../policy.html" title="Chapter&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">policy
  71         documentation for more details</a>.
  72       </p>
  73 <div class="blockquote"><blockquote class="blockquote"><p>
  74           <span class="inlinemediaobject"><img src="../../../graphs/lgamma.svg" align="middle"></span>
  75
  76         </p></blockquote></div>
  77 <p>
  78         The return type of these functions is computed using the <a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>result
  79         type calculation rules</em></span></a>: the result is of type <code class="computeroutput"><span class="keyword">double</span></code> if T is an integer type, or type T
  80         otherwise.
  81       </p>
  82 <h5>
  83 <a name="math_toolkit.sf_gamma.lgamma.h2"></a>
  84         <span class="phrase"><a name="math_toolkit.sf_gamma.lgamma.accuracy"></a></span><a class="link" href="lgamma.html#math_toolkit.sf_gamma.lgamma.accuracy">Accuracy</a>
  85       </h5>
  86 <p>
  87         The following table shows the peak errors (in units of epsilon) found on
  88         various platforms with various floating point types, along with comparisons
  89         to various other libraries. Unless otherwise specified any floating point
  90         type that is narrower than the one shown will have <a class="link" href="../relative_error.html#math_toolkit.relative_error.zero_error">effectively
  91         zero error</a>.
  92       </p>
  93 <p>
  94         Note that while the relative errors near the positive roots of lgamma are
  95         very low, the lgamma function has an infinite number of irrational roots
  96         for negative arguments: very close to these negative roots only a low absolute
  97         error can be guaranteed.
  98       </p>
  99 <div class="table">
 100 <a name="math_toolkit.sf_gamma.lgamma.table_lgamma"></a><p class="title"><b>Table&#160;8.3.&#160;Error rates for lgamma</b></p>
 101 <div class="table-contents"><table class="table" summary="Error rates for lgamma">
 102 <colgroup>
 103 <col>
 104 <col>
 105 <col>
 106 <col>
 107 <col>
 108 </colgroup>
 109 <thead><tr>
 110 <th>
 111               </th>
 112 <th>
 113                 <p>
 114                   GNU C++ version 7.1.0<br> linux<br> double
 115                 </p>
 116               </th>
 117 <th>
 118                 <p>
 119                   GNU C++ version 7.1.0<br> linux<br> long double
 120                 </p>
 121               </th>
 122 <th>
 123                 <p>
 124                   Sun compiler version 0x5150<br> Sun Solaris<br> long double
 125                 </p>
 126               </th>
 127 <th>
 128                 <p>
 129                   Microsoft Visual C++ version 14.1<br> Win32<br> double
 130                 </p>
 131               </th>
 132 </tr></thead>
 133 <tbody>
 134 <tr>
 135 <td>
 136                 <p>
 137                   factorials
 138                 </p>
 139               </td>
 140 <td>
 141                 <p>
 142                   <span class="blue">Max = 0&#949; (Mean = 0&#949;)</span><br> <br> (<span class="emphasis"><em>GSL
 143                   2.1:</em></span> Max = 33.6&#949; (Mean = 2.78&#949;))<br> (<span class="emphasis"><em>Rmath
 144                   3.2.3:</em></span> Max = 1.55&#949; (Mean = 0.592&#949;))
 145                 </p>
 146               </td>
 147 <td>
 148                 <p>
 149                   <span class="blue">Max = 0.991&#949; (Mean = 0.308&#949;)</span><br> <br>
 150                   (<span class="emphasis"><em>&lt;cmath&gt;:</em></span> Max = 1.67&#949; (Mean = 0.487&#949;))<br>
 151                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 1.67&#949; (Mean = 0.487&#949;))
 152                 </p>
 153               </td>
 154 <td>
 155                 <p>
 156                   <span class="blue">Max = 0.991&#949; (Mean = 0.383&#949;)</span><br> <br>
 157                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 1.36&#949; (Mean = 0.476&#949;))
 158                 </p>
 159               </td>
 160 <td>
 161                 <p>
 162                   <span class="blue">Max = 0.914&#949; (Mean = 0.175&#949;)</span><br> <br>
 163                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.958&#949; (Mean = 0.38&#949;))
 164                 </p>
 165               </td>
 166 </tr>
 167 <tr>
 168 <td>
 169                 <p>
 170                   near 0
 171                 </p>
 172               </td>
 173 <td>
 174                 <p>
 175                   <span class="blue">Max = 0&#949; (Mean = 0&#949;)</span><br> <br> (<span class="emphasis"><em>GSL
 176                   2.1:</em></span> Max = 5.21&#949; (Mean = 1.57&#949;))<br> (<span class="emphasis"><em>Rmath
 177                   3.2.3:</em></span> Max = 0&#949; (Mean = 0&#949;))
 178                 </p>
 179               </td>
 180 <td>
 181                 <p>
 182                   <span class="blue">Max = 1.42&#949; (Mean = 0.566&#949;)</span><br> <br>
 183                   (<span class="emphasis"><em>&lt;cmath&gt;:</em></span> Max = 0.964&#949; (Mean = 0.543&#949;))<br>
 184                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.964&#949; (Mean = 0.543&#949;))
 185                 </p>
 186               </td>
 187 <td>
 188                 <p>
 189                   <span class="blue">Max = 1.42&#949; (Mean = 0.566&#949;)</span><br> <br>
 190                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.964&#949; (Mean = 0.543&#949;))
 191                 </p>
 192               </td>
 193 <td>
 194                 <p>
 195                   <span class="blue">Max = 0.964&#949; (Mean = 0.462&#949;)</span><br> <br>
 196                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.962&#949; (Mean = 0.372&#949;))
 197                 </p>
 198               </td>
 199 </tr>
 200 <tr>
 201 <td>
 202                 <p>
 203                   near 1
 204                 </p>
 205               </td>
 206 <td>
 207                 <p>
 208                   <span class="blue">Max = 0&#949; (Mean = 0&#949;)</span><br> <br> (<span class="emphasis"><em>GSL
 209                   2.1:</em></span> Max = 442&#949; (Mean = 88.8&#949;))<br> (<span class="emphasis"><em>Rmath
 210                   3.2.3:</em></span> Max = 7.99e+04&#949; (Mean = 1.68e+04&#949;))
 211                 </p>
 212               </td>
 213 <td>
 214                 <p>
 215                   <span class="blue">Max = 0.948&#949; (Mean = 0.36&#949;)</span><br> <br>
 216                   (<span class="emphasis"><em>&lt;cmath&gt;:</em></span> Max = 0.615&#949; (Mean = 0.096&#949;))<br>
 217                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.615&#949; (Mean = 0.096&#949;))
 218                 </p>
 219               </td>
 220 <td>
 221                 <p>
 222                   <span class="blue">Max = 0.948&#949; (Mean = 0.36&#949;)</span><br> <br>
 223                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 1.71&#949; (Mean = 0.581&#949;))
 224                 </p>
 225               </td>
 226 <td>
 227                 <p>
 228                   <span class="blue">Max = 0.867&#949; (Mean = 0.468&#949;)</span><br> <br>
 229                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.906&#949; (Mean = 0.565&#949;))
 230                 </p>
 231               </td>
 232 </tr>
 233 <tr>
 234 <td>
 235                 <p>
 236                   near 2
 237                 </p>
 238               </td>
 239 <td>
 240                 <p>
 241                   <span class="blue">Max = 0&#949; (Mean = 0&#949;)</span><br> <br> (<span class="emphasis"><em>GSL
 242                   2.1:</em></span> Max = 1.17e+03&#949; (Mean = 274&#949;))<br> (<span class="emphasis"><em>Rmath
 243                   3.2.3:</em></span> Max = 2.63e+05&#949; (Mean = 5.84e+04&#949;))
 244                 </p>
 245               </td>
 246 <td>
 247                 <p>
 248                   <span class="blue">Max = 0.878&#949; (Mean = 0.242&#949;)</span><br> <br>
 249                   (<span class="emphasis"><em>&lt;cmath&gt;:</em></span> Max = 0.741&#949; (Mean = 0.263&#949;))<br>
 250                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.741&#949; (Mean = 0.263&#949;))
 251                 </p>
 252               </td>
 253 <td>
 254                 <p>
 255                   <span class="blue">Max = 0.878&#949; (Mean = 0.242&#949;)</span><br> <br>
 256                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.598&#949; (Mean = 0.235&#949;))
 257                 </p>
 258               </td>
 259 <td>
 260                 <p>
 261                   <span class="blue">Max = 0.591&#949; (Mean = 0.159&#949;)</span><br> <br>
 262                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.741&#949; (Mean = 0.473&#949;))
 263                 </p>
 264               </td>
 265 </tr>
 266 <tr>
 267 <td>
 268                 <p>
 269                   near -10
 270                 </p>
 271               </td>
 272 <td>
 273                 <p>
 274                   <span class="blue">Max = 0&#949; (Mean = 0&#949;)</span><br> <br> (<span class="emphasis"><em>GSL
 275                   2.1:</em></span> Max = 24.9&#949; (Mean = 4.6&#949;))<br> (<span class="emphasis"><em>Rmath
 276                   3.2.3:</em></span> Max = 4.22&#949; (Mean = 1.26&#949;))
 277                 </p>
 278               </td>
 279 <td>
 280                 <p>
 281                   <span class="blue">Max = 3.81&#949; (Mean = 1.01&#949;)</span><br> <br>
 282                   (<span class="emphasis"><em>&lt;cmath&gt;:</em></span> Max = 0.997&#949; (Mean = 0.412&#949;))<br>
 283                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.997&#949; (Mean = 0.412&#949;))
 284                 </p>
 285               </td>
 286 <td>
 287                 <p>
 288                   <span class="blue">Max = 3.81&#949; (Mean = 1.01&#949;)</span><br> <br>
 289                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 3.04&#949; (Mean = 1.01&#949;))
 290                 </p>
 291               </td>
 292 <td>
 293                 <p>
 294                   <span class="blue">Max = 4.22&#949; (Mean = 1.33&#949;)</span><br> <br>
 295                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.997&#949; (Mean = 0.444&#949;))
 296                 </p>
 297               </td>
 298 </tr>
 299 <tr>
 300 <td>
 301                 <p>
 302                   near -55
 303                 </p>
 304               </td>
 305 <td>
 306                 <p>
 307                   <span class="blue">Max = 0&#949; (Mean = 0&#949;)</span><br> <br> (<span class="emphasis"><em>GSL
 308                   2.1:</em></span> Max = 7.02&#949; (Mean = 1.47&#949;))<br> (<span class="emphasis"><em>Rmath
 309                   3.2.3:</em></span> Max = 250&#949; (Mean = 60.9&#949;))
 310                 </p>
 311               </td>
 312 <td>
 313                 <p>
 314                   <span class="blue">Max = 0.821&#949; (Mean = 0.513&#949;)</span><br> <br>
 315                   (<span class="emphasis"><em>&lt;cmath&gt;:</em></span> Max = 1.58&#949; (Mean = 0.672&#949;))<br>
 316                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 1.58&#949; (Mean = 0.672&#949;))
 317                 </p>
 318               </td>
 319 <td>
 320                 <p>
 321                   <span class="blue">Max = 1.59&#949; (Mean = 0.587&#949;)</span><br> <br>
 322                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 0.821&#949; (Mean = 0.674&#949;))
 323                 </p>
 324               </td>
 325 <td>
 326                 <p>
 327                   <span class="blue">Max = 0.821&#949; (Mean = 0.419&#949;)</span><br> <br>
 328                   (<span class="emphasis"><em>&lt;math.h&gt;:</em></span> Max = 249&#949; (Mean = 43.1&#949;))
 329                 </p>
 330               </td>
 331 </tr>
 332 </tbody>
 333 </table></div>
 334 </div>
 335 <br class="table-break"><p>
 336         The following error plot are based on an exhaustive search of the functions
 337         domain, MSVC-15.5 at <code class="computeroutput"><span class="keyword">double</span></code>
 338         precision, and GCC-7.1/Ubuntu for <code class="computeroutput"><span class="keyword">long</span>
 339         <span class="keyword">double</span></code> and <code class="computeroutput"><span class="identifier">__float128</span></code>.
 340       </p>
 341 <div class="blockquote"><blockquote class="blockquote"><p>
 342           <span class="inlinemediaobject"><img src="../../../graphs/lgamma__double.svg" align="middle"></span>
 343
 344         </p></blockquote></div>
 345 <div class="blockquote"><blockquote class="blockquote"><p>
 346           <span class="inlinemediaobject"><img src="../../../graphs/lgamma__80_bit_long_double.svg" align="middle"></span>
 347
 348         </p></blockquote></div>
 349 <div class="blockquote"><blockquote class="blockquote"><p>
 350           <span class="inlinemediaobject"><img src="../../../graphs/lgamma____float128.svg" align="middle"></span>
 351
 352         </p></blockquote></div>
 353 <h5>
 354 <a name="math_toolkit.sf_gamma.lgamma.h3"></a>
 355         <span class="phrase"><a name="math_toolkit.sf_gamma.lgamma.testing"></a></span><a class="link" href="lgamma.html#math_toolkit.sf_gamma.lgamma.testing">Testing</a>
 356       </h5>
 357 <p>
 358         The main tests for this function involve comparisons against the logs of
 359         the factorials which can be independently calculated to very high accuracy.
 360       </p>
 361 <p>
 362         Random tests in key problem areas are also used.
 363       </p>
 364 <h5>
 365 <a name="math_toolkit.sf_gamma.lgamma.h4"></a>
 366         <span class="phrase"><a name="math_toolkit.sf_gamma.lgamma.implementation"></a></span><a class="link" href="lgamma.html#math_toolkit.sf_gamma.lgamma.implementation">Implementation</a>
 367       </h5>
 368 <p>
 369         The generic version of this function is implemented using Sterling's approximation
 370         for large arguments:
 371       </p>
 372 <div class="blockquote"><blockquote class="blockquote"><p>
 373           <span class="inlinemediaobject"><img src="../../../equations/gamma6.svg"></span>
 374
 375         </p></blockquote></div>
 376 <p>
 377         For small arguments, the logarithm of tgamma is used.
 378       </p>
 379 <p>
 380         For negative <span class="emphasis"><em>z</em></span> the logarithm version of the reflection
 381         formula is used:
 382       </p>
 383 <div class="blockquote"><blockquote class="blockquote"><p>
 384           <span class="inlinemediaobject"><img src="../../../equations/lgamm3.svg"></span>
 385
 386         </p></blockquote></div>
 387 <p>
 388         For types of known precision, the <a class="link" href="../lanczos.html" title="The Lanczos Approximation">Lanczos
 389         approximation</a> is used, a traits class <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">lanczos</span><span class="special">::</span><span class="identifier">lanczos_traits</span></code>
 390         maps type T to an appropriate approximation. The logarithmic version of the
 391         <a class="link" href="../lanczos.html" title="The Lanczos Approximation">Lanczos approximation</a> is:
 392       </p>
 393 <div class="blockquote"><blockquote class="blockquote"><p>
 394           <span class="inlinemediaobject"><img src="../../../equations/lgamm4.svg"></span>
 395
 396         </p></blockquote></div>
 397 <p>
 398         Where L<sub>e,g</sub> is the Lanczos sum, scaled by e<sup>g</sup>.
 399       </p>
 400 <p>
 401         As before the reflection formula is used for <span class="emphasis"><em>z &lt; 0</em></span>.
 402       </p>
 403 <p>
 404         When z is very near 1 or 2, then the logarithmic version of the <a class="link" href="../lanczos.html" title="The Lanczos Approximation">Lanczos
 405         approximation</a> suffers very badly from cancellation error: indeed for
 406         values sufficiently close to 1 or 2, arbitrarily large relative errors can
 407         be obtained (even though the absolute error is tiny).
 408       </p>
 409 <p>
 410         For types with up to 113 bits of precision (up to and including 128-bit long
 411         doubles), root-preserving rational approximations <a class="link" href="../sf_implementation.html#math_toolkit.sf_implementation.rational_approximations_used">devised
 412         by JM</a> are used over the intervals [1,2] and [2,3]. Over the interval
 413         [2,3] the approximation form used is:
 414       </p>
 415 <pre class="programlisting"><span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">z</span><span class="special">-</span><span class="number">2</span><span class="special">)(</span><span class="identifier">z</span><span class="special">+</span><span class="number">1</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">z</span><span class="special">-</span><span class="number">2</span><span class="special">));</span>
 416 </pre>
 417 <p>
 418         Where Y is a constant, and R(z-2) is the rational approximation: optimised
 419         so that its absolute error is tiny compared to Y. In addition, small values
 420         of z greater than 3 can handled by argument reduction using the recurrence
 421         relation:
 422       </p>
 423 <pre class="programlisting"><span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">z</span><span class="special">+</span><span class="number">1</span><span class="special">)</span> <span class="special">=</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special">+</span> <span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span>
 424 </pre>
 425 <p>
 426         Over the interval [1,2] two approximations have to be used, one for small
 427         z uses:
 428       </p>
 429 <pre class="programlisting"><span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">z</span><span class="special">-</span><span class="number">1</span><span class="special">)(</span><span class="identifier">z</span><span class="special">-</span><span class="number">2</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">z</span><span class="special">-</span><span class="number">1</span><span class="special">));</span>
 430 </pre>
 431 <p>
 432         Once again Y is a constant, and R(z-1) is optimised for low absolute error
 433         compared to Y. For z &gt; 1.5 the above form wouldn't converge to a minimax
 434         solution but this similar form does:
 435       </p>
 436 <pre class="programlisting"><span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special">=</span> <span class="special">(</span><span class="number">2</span><span class="special">-</span><span class="identifier">z</span><span class="special">)(</span><span class="number">1</span><span class="special">-</span><span class="identifier">z</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="number">2</span><span class="special">-</span><span class="identifier">z</span><span class="special">));</span>
 437 </pre>
 438 <p>
 439         Finally for z &lt; 1 the recurrence relation can be used to move to z &gt;
 440         1:
 441       </p>
 442 <pre class="programlisting"><span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special">=</span> <span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">z</span><span class="special">+</span><span class="number">1</span><span class="special">)</span> <span class="special">-</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span>
 443 </pre>
 444 <p>
 445         Note that while this involves a subtraction, it appears not to suffer from
 446         cancellation error: as z decreases from 1 the <code class="computeroutput"><span class="special">-</span><span class="identifier">log</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span></code> term grows positive much more rapidly than
 447         the <code class="computeroutput"><span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">z</span><span class="special">+</span><span class="number">1</span><span class="special">)</span></code> term becomes negative. So in this specific
 448         case, significant digits are preserved, rather than cancelled.
 449       </p>
 450 <p>
 451         For other types which do have a <a class="link" href="../lanczos.html" title="The Lanczos Approximation">Lanczos
 452         approximation</a> defined for them the current solution is as follows:
 453         imagine we balance the two terms in the <a class="link" href="../lanczos.html" title="The Lanczos Approximation">Lanczos
 454         approximation</a> by dividing the power term by its value at <span class="emphasis"><em>z
 455         = 1</em></span>, and then multiplying the Lanczos coefficients by the same
 456         value. Now each term will take the value 1 at <span class="emphasis"><em>z = 1</em></span>
 457         and we can rearrange the power terms in terms of log1p. Likewise if we subtract
 458         1 from the Lanczos sum part (algebraically, by subtracting the value of each
 459         term at <span class="emphasis"><em>z = 1</em></span>), we obtain a new summation that can be
 460         also be fed into log1p. Crucially, all of the terms tend to zero, as <span class="emphasis"><em>z
 461         -&gt; 1</em></span>:
 462       </p>
 463 <div class="blockquote"><blockquote class="blockquote"><p>
 464           <span class="inlinemediaobject"><img src="../../../equations/lgamm5.svg"></span>
 465
 466         </p></blockquote></div>
 467 <p>
 468         The C<sub>k</sub> terms in the above are the same as in the <a class="link" href="../lanczos.html" title="The Lanczos Approximation">Lanczos
 469         approximation</a>.
 470       </p>
 471 <p>
 472         A similar rearrangement can be performed at <span class="emphasis"><em>z = 2</em></span>:
 473       </p>
 474 <div class="blockquote"><blockquote class="blockquote"><p>
 475           <span class="inlinemediaobject"><img src="../../../equations/lgamm6.svg"></span>
 476
 477         </p></blockquote></div>
 478 </div>
 479 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 480 <td align="left"></td>
 481 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2019 Nikhar
 482       Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
 483       Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
 484       R&#229;de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
 485       Daryle Walker and Xiaogang Zhang<p>
 486         Distributed under the Boost Software License, Version 1.0. (See accompanying
 487         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>)
 488       </p>
 489 </div></td>
 490 </tr></table>
 491 <hr>
 492 <div class="spirit-nav">
 493 <a accesskey="p" href="tgamma.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="digamma.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
 494 </div>
 495 </body>
 496 </html>