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