1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 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.12.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"><</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">></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">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"><</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">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 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> 42 43<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></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"><</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> 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 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</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 21. 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 21. 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 8.3. 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ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 143 2.1:</em></span> Max = 33.6ε (Mean = 2.78ε))<br> (<span class="emphasis"><em>Rmath 144 3.2.3:</em></span> Max = 1.55ε (Mean = 0.592ε)) 145 </p> 146 </td> 147<td> 148 <p> 149 <span class="blue">Max = 0.991ε (Mean = 0.308ε)</span><br> <br> 150 (<span class="emphasis"><em><cmath>:</em></span> Max = 1.67ε (Mean = 0.487ε))<br> 151 (<span class="emphasis"><em><math.h>:</em></span> Max = 1.67ε (Mean = 0.487ε)) 152 </p> 153 </td> 154<td> 155 <p> 156 <span class="blue">Max = 0.991ε (Mean = 0.383ε)</span><br> <br> 157 (<span class="emphasis"><em><math.h>:</em></span> Max = 1.36ε (Mean = 0.476ε)) 158 </p> 159 </td> 160<td> 161 <p> 162 <span class="blue">Max = 0.914ε (Mean = 0.175ε)</span><br> <br> 163 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.958ε (Mean = 0.38ε)) 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ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 176 2.1:</em></span> Max = 5.21ε (Mean = 1.57ε))<br> (<span class="emphasis"><em>Rmath 177 3.2.3:</em></span> Max = 0ε (Mean = 0ε)) 178 </p> 179 </td> 180<td> 181 <p> 182 <span class="blue">Max = 1.42ε (Mean = 0.566ε)</span><br> <br> 183 (<span class="emphasis"><em><cmath>:</em></span> Max = 0.964ε (Mean = 0.543ε))<br> 184 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.964ε (Mean = 0.543ε)) 185 </p> 186 </td> 187<td> 188 <p> 189 <span class="blue">Max = 1.42ε (Mean = 0.566ε)</span><br> <br> 190 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.964ε (Mean = 0.543ε)) 191 </p> 192 </td> 193<td> 194 <p> 195 <span class="blue">Max = 0.964ε (Mean = 0.462ε)</span><br> <br> 196 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.962ε (Mean = 0.372ε)) 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ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 209 2.1:</em></span> Max = 442ε (Mean = 88.8ε))<br> (<span class="emphasis"><em>Rmath 210 3.2.3:</em></span> Max = 7.99e+04ε (Mean = 1.68e+04ε)) 211 </p> 212 </td> 213<td> 214 <p> 215 <span class="blue">Max = 0.948ε (Mean = 0.36ε)</span><br> <br> 216 (<span class="emphasis"><em><cmath>:</em></span> Max = 0.615ε (Mean = 0.096ε))<br> 217 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.615ε (Mean = 0.096ε)) 218 </p> 219 </td> 220<td> 221 <p> 222 <span class="blue">Max = 0.948ε (Mean = 0.36ε)</span><br> <br> 223 (<span class="emphasis"><em><math.h>:</em></span> Max = 1.71ε (Mean = 0.581ε)) 224 </p> 225 </td> 226<td> 227 <p> 228 <span class="blue">Max = 0.867ε (Mean = 0.468ε)</span><br> <br> 229 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.906ε (Mean = 0.565ε)) 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ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 242 2.1:</em></span> Max = 1.17e+03ε (Mean = 274ε))<br> (<span class="emphasis"><em>Rmath 243 3.2.3:</em></span> Max = 2.63e+05ε (Mean = 5.84e+04ε)) 244 </p> 245 </td> 246<td> 247 <p> 248 <span class="blue">Max = 0.878ε (Mean = 0.242ε)</span><br> <br> 249 (<span class="emphasis"><em><cmath>:</em></span> Max = 0.741ε (Mean = 0.263ε))<br> 250 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.741ε (Mean = 0.263ε)) 251 </p> 252 </td> 253<td> 254 <p> 255 <span class="blue">Max = 0.878ε (Mean = 0.242ε)</span><br> <br> 256 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.598ε (Mean = 0.235ε)) 257 </p> 258 </td> 259<td> 260 <p> 261 <span class="blue">Max = 0.591ε (Mean = 0.159ε)</span><br> <br> 262 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.741ε (Mean = 0.473ε)) 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ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 275 2.1:</em></span> Max = 24.9ε (Mean = 4.6ε))<br> (<span class="emphasis"><em>Rmath 276 3.2.3:</em></span> Max = 4.22ε (Mean = 1.26ε)) 277 </p> 278 </td> 279<td> 280 <p> 281 <span class="blue">Max = 3.81ε (Mean = 1.01ε)</span><br> <br> 282 (<span class="emphasis"><em><cmath>:</em></span> Max = 0.997ε (Mean = 0.412ε))<br> 283 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.997ε (Mean = 0.412ε)) 284 </p> 285 </td> 286<td> 287 <p> 288 <span class="blue">Max = 3.81ε (Mean = 1.01ε)</span><br> <br> 289 (<span class="emphasis"><em><math.h>:</em></span> Max = 3.04ε (Mean = 1.01ε)) 290 </p> 291 </td> 292<td> 293 <p> 294 <span class="blue">Max = 4.22ε (Mean = 1.33ε)</span><br> <br> 295 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.997ε (Mean = 0.444ε)) 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ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 308 2.1:</em></span> Max = 7.02ε (Mean = 1.47ε))<br> (<span class="emphasis"><em>Rmath 309 3.2.3:</em></span> Max = 250ε (Mean = 60.9ε)) 310 </p> 311 </td> 312<td> 313 <p> 314 <span class="blue">Max = 0.821ε (Mean = 0.513ε)</span><br> <br> 315 (<span class="emphasis"><em><cmath>:</em></span> Max = 1.58ε (Mean = 0.672ε))<br> 316 (<span class="emphasis"><em><math.h>:</em></span> Max = 1.58ε (Mean = 0.672ε)) 317 </p> 318 </td> 319<td> 320 <p> 321 <span class="blue">Max = 1.59ε (Mean = 0.587ε)</span><br> <br> 322 (<span class="emphasis"><em><math.h>:</em></span> Max = 0.821ε (Mean = 0.674ε)) 323 </p> 324 </td> 325<td> 326 <p> 327 <span class="blue">Max = 0.821ε (Mean = 0.419ε)</span><br> <br> 328 (<span class="emphasis"><em><math.h>:</em></span> Max = 249ε (Mean = 43.1ε)) 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 < 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 > 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 < 1 the recurrence relation can be used to move to z > 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 -> 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 © 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å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> 497
