| <html> |
| <head> |
| <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII"> |
| <title>The Lanczos Approximation</title> |
| <link rel="stylesheet" href="../math.css" type="text/css"> |
| <meta name="generator" content="DocBook XSL Stylesheets V1.77.1"> |
| <link rel="home" href="../index.html" title="Math Toolkit 2.2.0"> |
| <link rel="up" href="../backgrounders.html" title="Chapter 16. Backgrounders"> |
| <link rel="prev" href="relative_error.html" title="Relative Error"> |
| <link rel="next" href="remez.html" title="The Remez Method"> |
| </head> |
| <body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> |
| <table cellpadding="2" width="100%"><tr> |
| <td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td> |
| <td align="center"><a href="../../../../../index.html">Home</a></td> |
| <td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td> |
| <td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> |
| <td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> |
| <td align="center"><a href="../../../../../more/index.htm">More</a></td> |
| </tr></table> |
| <hr> |
| <div class="spirit-nav"> |
| <a accesskey="p" href="relative_error.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../backgrounders.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="remez.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> |
| </div> |
| <div class="section"> |
| <div class="titlepage"><div><div><h2 class="title" style="clear: both"> |
| <a name="math_toolkit.lanczos"></a><a class="link" href="lanczos.html" title="The Lanczos Approximation">The Lanczos Approximation</a> |
| </h2></div></div></div> |
| <h5> |
| <a name="math_toolkit.lanczos.h0"></a> |
| <span class="phrase"><a name="math_toolkit.lanczos.motivation"></a></span><a class="link" href="lanczos.html#math_toolkit.lanczos.motivation">Motivation</a> |
| </h5> |
| <p> |
| <span class="emphasis"><em>Why base gamma and gamma-like functions on the Lanczos approximation?</em></span> |
| </p> |
| <p> |
| First of all I should make clear that for the gamma function over real numbers |
| (as opposed to complex ones) the Lanczos approximation (See <a href="http://en.wikipedia.org/wiki/Lanczos_approximation" target="_top">Wikipedia |
| or </a> <a href="http://mathworld.wolfram.com/LanczosApproximation.html" target="_top">Mathworld</a>) |
| appears to offer no clear advantage over more traditional methods such as |
| <a href="http://en.wikipedia.org/wiki/Stirling_approximation" target="_top">Stirling's |
| approximation</a>. <a class="link" href="lanczos.html#pugh">Pugh</a> carried out an extensive |
| comparison of the various methods available and discovered that they were all |
| very similar in terms of complexity and relative error. However, the Lanczos |
| approximation does have a couple of properties that make it worthy of further |
| consideration: |
| </p> |
| <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> |
| <li class="listitem"> |
| The approximation has an easy to compute truncation error that holds for |
| all <span class="emphasis"><em>z > 0</em></span>. In practice that means we can use the |
| same approximation for all <span class="emphasis"><em>z > 0</em></span>, and be certain |
| that no matter how large or small <span class="emphasis"><em>z</em></span> is, the truncation |
| error will <span class="emphasis"><em>at worst</em></span> be bounded by some finite value. |
| </li> |
| <li class="listitem"> |
| The approximation has a form that is particularly amenable to analytic |
| manipulation, in particular ratios of gamma or gamma-like functions are |
| particularly easy to compute without resorting to logarithms. |
| </li> |
| </ul></div> |
| <p> |
| It is the combination of these two properties that make the approximation attractive: |
| Stirling's approximation is highly accurate for large z, and has some of the |
| same analytic properties as the Lanczos approximation, but can't easily be |
| used across the whole range of z. |
| </p> |
| <p> |
| As the simplest example, consider the ratio of two gamma functions: one could |
| compute the result via lgamma: |
| </p> |
| <pre class="programlisting"><span class="identifier">exp</span><span class="special">(</span><span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">a</span><span class="special">)</span> <span class="special">-</span> <span class="identifier">lgamma</span><span class="special">(</span><span class="identifier">b</span><span class="special">));</span> |
| </pre> |
| <p> |
| However, even if lgamma is uniformly accurate to 0.5ulp, the worst case relative |
| error in the above can easily be shown to be: |
| </p> |
| <pre class="programlisting"><span class="identifier">Erel</span> <span class="special">></span> <span class="identifier">a</span> <span class="special">*</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">a</span><span class="special">)/</span><span class="number">2</span> <span class="special">+</span> <span class="identifier">b</span> <span class="special">*</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">b</span><span class="special">)/</span><span class="number">2</span> |
| </pre> |
| <p> |
| For small <span class="emphasis"><em>a</em></span> and <span class="emphasis"><em>b</em></span> that's not a problem, |
| but to put the relationship another way: <span class="emphasis"><em>each time a and b increase |
| in magnitude by a factor of 10, at least one decimal digit of precision will |
| be lost.</em></span> |
| </p> |
| <p> |
| In contrast, by analytically combining like power terms in a ratio of Lanczos |
| approximation's, these errors can be virtually eliminated for small <span class="emphasis"><em>a</em></span> |
| and <span class="emphasis"><em>b</em></span>, and kept under control for very large (or very |
| small for that matter) <span class="emphasis"><em>a</em></span> and <span class="emphasis"><em>b</em></span>. Of |
| course, computing large powers is itself a notoriously hard problem, but even |
| so, analytic combinations of Lanczos approximations can make the difference |
| between obtaining a valid result, or simply garbage. Refer to the implementation |
| notes for the <a class="link" href="sf_beta/beta_function.html" title="Beta">beta</a> |
| function for an example of this method in practice. The incomplete <a class="link" href="sf_gamma/igamma.html" title="Incomplete Gamma Functions">gamma_p |
| gamma</a> and <a class="link" href="sf_beta/ibeta_function.html" title="Incomplete Beta Functions">beta</a> |
| functions use similar analytic combinations of power terms, to combine gamma |
| and beta functions divided by large powers into single (simpler) expressions. |
| </p> |
| <h5> |
| <a name="math_toolkit.lanczos.h1"></a> |
| <span class="phrase"><a name="math_toolkit.lanczos.the_approximation"></a></span><a class="link" href="lanczos.html#math_toolkit.lanczos.the_approximation">The |
| Approximation</a> |
| </h5> |
| <p> |
| The Lanczos Approximation to the Gamma Function is given by: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos0.svg"></span> |
| </p> |
| <p> |
| Where S<sub>g</sub>(z) is an infinite sum, that is convergent for all z > 0, and <span class="emphasis"><em>g</em></span> |
| is an arbitrary parameter that controls the "shape" of the terms |
| in the sum which is given by: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos0a.svg"></span> |
| </p> |
| <p> |
| With individual coefficients defined in closed form by: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos0b.svg"></span> |
| </p> |
| <p> |
| However, evaluation of the sum in that form can lead to numerical instability |
| in the computation of the ratios of rising and falling factorials (effectively |
| we're multiplying by a series of numbers very close to 1, so roundoff errors |
| can accumulate quite rapidly). |
| </p> |
| <p> |
| The Lanczos approximation is therefore often written in partial fraction form |
| with the leading constants absorbed by the coefficients in the sum: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos1.svg"></span> |
| </p> |
| <p> |
| where: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos2.svg"></span> |
| </p> |
| <p> |
| Again parameter <span class="emphasis"><em>g</em></span> is an arbitrarily chosen constant, and |
| <span class="emphasis"><em>N</em></span> is an arbitrarily chosen number of terms to evaluate |
| in the "Lanczos sum" part. |
| </p> |
| <div class="note"><table border="0" summary="Note"> |
| <tr> |
| <td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../doc/src/images/note.png"></td> |
| <th align="left">Note</th> |
| </tr> |
| <tr><td align="left" valign="top"><p> |
| Some authors choose to define the sum from k=1 to N, and hence end up with |
| N+1 coefficients. This happens to confuse both the following discussion and |
| the code (since C++ deals with half open array ranges, rather than the closed |
| range of the sum). This convention is consistent with <a class="link" href="lanczos.html#godfrey">Godfrey</a>, |
| but not <a class="link" href="lanczos.html#pugh">Pugh</a>, so take care when referring to |
| the literature in this field. |
| </p></td></tr> |
| </table></div> |
| <h5> |
| <a name="math_toolkit.lanczos.h2"></a> |
| <span class="phrase"><a name="math_toolkit.lanczos.computing_the_coefficients"></a></span><a class="link" href="lanczos.html#math_toolkit.lanczos.computing_the_coefficients">Computing |
| the Coefficients</a> |
| </h5> |
| <p> |
| The coefficients C0..CN-1 need to be computed from <span class="emphasis"><em>N</em></span> and |
| <span class="emphasis"><em>g</em></span> at high precision, and then stored as part of the program. |
| Calculation of the coefficients is performed via the method of <a class="link" href="lanczos.html#godfrey">Godfrey</a>; |
| let the constants be contained in a column vector P, then: |
| </p> |
| <p> |
| P = D B C F |
| </p> |
| <p> |
| where B is an NxN matrix: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos4.svg"></span> |
| </p> |
| <p> |
| D is an NxN matrix: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos3.svg"></span> |
| </p> |
| <p> |
| C is an NxN matrix: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos5.svg"></span> |
| </p> |
| <p> |
| and F is an N element column vector: |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos6.svg"></span> |
| </p> |
| <p> |
| Note than the matrices B, D and C contain all integer terms and depend only |
| on <span class="emphasis"><em>N</em></span>, this product should be computed first, and then |
| multiplied by <span class="emphasis"><em>F</em></span> as the last step. |
| </p> |
| <h5> |
| <a name="math_toolkit.lanczos.h3"></a> |
| <span class="phrase"><a name="math_toolkit.lanczos.choosing_the_right_parameters"></a></span><a class="link" href="lanczos.html#math_toolkit.lanczos.choosing_the_right_parameters">Choosing |
| the Right Parameters</a> |
| </h5> |
| <p> |
| The trick is to choose <span class="emphasis"><em>N</em></span> and <span class="emphasis"><em>g</em></span> to |
| give the desired level of accuracy: choosing a small value for <span class="emphasis"><em>g</em></span> |
| leads to a strictly convergent series, but one which converges only slowly. |
| Choosing a larger value of <span class="emphasis"><em>g</em></span> causes the terms in the series |
| to be large and/or divergent for about the first <span class="emphasis"><em>g-1</em></span> terms, |
| and to then suddenly converge with a "crunch". |
| </p> |
| <p> |
| <a class="link" href="lanczos.html#pugh">Pugh</a> has determined the optimal value of <span class="emphasis"><em>g</em></span> |
| for <span class="emphasis"><em>N</em></span> in the range <span class="emphasis"><em>1 <= N <= 60</em></span>: |
| unfortunately in practice choosing these values leads to cancellation errors |
| in the Lanczos sum as the largest term in the (alternating) series is approximately |
| 1000 times larger than the result. These optimal values appear not to be useful |
| in practice unless the evaluation can be done with a number of guard digits |
| <span class="emphasis"><em>and</em></span> the coefficients are stored at higher precision than |
| that desired in the result. These values are best reserved for say, computing |
| to float precision with double precision arithmetic. |
| </p> |
| <div class="table"> |
| <a name="math_toolkit.lanczos.optimal_choices_for_n_and_g_when"></a><p class="title"><b>Table 16.1. Optimal choices for N and g when computing with guard digits (source: |
| Pugh)</b></p> |
| <div class="table-contents"><table class="table" summary="Optimal choices for N and g when computing with guard digits (source: |
| Pugh)"> |
| <colgroup> |
| <col> |
| <col> |
| <col> |
| <col> |
| </colgroup> |
| <thead><tr> |
| <th> |
| <p> |
| Significand Size |
| </p> |
| </th> |
| <th> |
| <p> |
| N |
| </p> |
| </th> |
| <th> |
| <p> |
| g |
| </p> |
| </th> |
| <th> |
| <p> |
| Max Error |
| </p> |
| </th> |
| </tr></thead> |
| <tbody> |
| <tr> |
| <td> |
| <p> |
| 24 |
| </p> |
| </td> |
| <td> |
| <p> |
| 6 |
| </p> |
| </td> |
| <td> |
| <p> |
| 5.581 |
| </p> |
| </td> |
| <td> |
| <p> |
| 9.51e-12 |
| </p> |
| </td> |
| </tr> |
| <tr> |
| <td> |
| <p> |
| 53 |
| </p> |
| </td> |
| <td> |
| <p> |
| 13 |
| </p> |
| </td> |
| <td> |
| <p> |
| 13.144565 |
| </p> |
| </td> |
| <td> |
| <p> |
| 9.2213e-23 |
| </p> |
| </td> |
| </tr> |
| </tbody> |
| </table></div> |
| </div> |
| <br class="table-break"><p> |
| The alternative described by <a class="link" href="lanczos.html#godfrey">Godfrey</a> is to perform |
| an exhaustive search of the <span class="emphasis"><em>N</em></span> and <span class="emphasis"><em>g</em></span> |
| parameter space to determine the optimal combination for a given <span class="emphasis"><em>p</em></span> |
| digit floating-point type. Repeating this work found a good approximation for |
| double precision arithmetic (close to the one <a class="link" href="lanczos.html#godfrey">Godfrey</a> |
| found), but failed to find really good approximations for 80 or 128-bit long |
| doubles. Further it was observed that the approximations obtained tended to |
| optimised for the small values of z (1 < z < 200) used to test the implementation |
| against the factorials. Computing ratios of gamma functions with large arguments |
| were observed to suffer from error resulting from the truncation of the Lancozos |
| series. |
| </p> |
| <p> |
| <a class="link" href="lanczos.html#pugh">Pugh</a> identified all the locations where the theoretical |
| error of the approximation were at a minimum, but unfortunately has published |
| only the largest of these minima. However, he makes the observation that the |
| minima coincide closely with the location where the first neglected term (a<sub>N</sub>) |
| in the Lanczos series S<sub>g</sub>(z) changes sign. These locations are quite easy to |
| locate, albeit with considerable computer time. These "sweet spots" |
| need only be computed once, tabulated, and then searched when required for |
| an approximation that delivers the required precision for some fixed precision |
| type. |
| </p> |
| <p> |
| Unfortunately, following this path failed to find a really good approximation |
| for 128-bit long doubles, and those found for 64 and 80-bit reals required |
| an excessive number of terms. There are two competing issues here: high precision |
| requires a large value of <span class="emphasis"><em>g</em></span>, but avoiding cancellation |
| errors in the evaluation requires a small <span class="emphasis"><em>g</em></span>. |
| </p> |
| <p> |
| At this point note that the Lanczos sum can be converted into rational form |
| (a ratio of two polynomials, obtained from the partial-fraction form using |
| polynomial arithmetic), and doing so changes the coefficients so that <span class="emphasis"><em>they |
| are all positive</em></span>. That means that the sum in rational form can be |
| evaluated without cancellation error, albeit with double the number of coefficients |
| for a given N. Repeating the search of the "sweet spots", this time |
| evaluating the Lanczos sum in rational form, and testing only those "sweet |
| spots" whose theoretical error is less than the machine epsilon for the |
| type being tested, yielded good approximations for all the types tested. The |
| optimal values found were quite close to the best cases reported by <a class="link" href="lanczos.html#pugh">Pugh</a> |
| (just slightly larger <span class="emphasis"><em>N</em></span> and slightly smaller <span class="emphasis"><em>g</em></span> |
| for a given precision than <a class="link" href="lanczos.html#pugh">Pugh</a> reports), and even |
| though converting to rational form doubles the number of stored coefficients, |
| it should be noted that half of them are integers (and therefore require less |
| storage space) and the approximations require a smaller <span class="emphasis"><em>N</em></span> |
| than would otherwise be required, so fewer floating point operations may be |
| required overall. |
| </p> |
| <p> |
| The following table shows the optimal values for <span class="emphasis"><em>N</em></span> and |
| <span class="emphasis"><em>g</em></span> when computing at fixed precision. These should be taken |
| as work in progress: there are no values for 106-bit significand machines (Darwin |
| long doubles & NTL quad_float), and further optimisation of the values |
| of <span class="emphasis"><em>g</em></span> may be possible. Errors given in the table are estimates |
| of the error due to truncation of the Lanczos infinite series to <span class="emphasis"><em>N</em></span> |
| terms. They are calculated from the sum of the first five neglected terms - |
| and are known to be rather pessimistic estimates - although it is noticeable |
| that the best combinations of <span class="emphasis"><em>N</em></span> and <span class="emphasis"><em>g</em></span> |
| occurred when the estimated truncation error almost exactly matches the machine |
| epsilon for the type in question. |
| </p> |
| <div class="table"> |
| <a name="math_toolkit.lanczos.optimum_value_for_n_and_g_when_c"></a><p class="title"><b>Table 16.2. Optimum value for N and g when computing at fixed precision</b></p> |
| <div class="table-contents"><table class="table" summary="Optimum value for N and g when computing at fixed precision"> |
| <colgroup> |
| <col> |
| <col> |
| <col> |
| <col> |
| <col> |
| </colgroup> |
| <thead><tr> |
| <th> |
| <p> |
| Significand Size |
| </p> |
| </th> |
| <th> |
| <p> |
| Platform/Compiler Used |
| </p> |
| </th> |
| <th> |
| <p> |
| N |
| </p> |
| </th> |
| <th> |
| <p> |
| g |
| </p> |
| </th> |
| <th> |
| <p> |
| Max Truncation Error |
| </p> |
| </th> |
| </tr></thead> |
| <tbody> |
| <tr> |
| <td> |
| <p> |
| 24 |
| </p> |
| </td> |
| <td> |
| <p> |
| Win32, VC++ 7.1 |
| </p> |
| </td> |
| <td> |
| <p> |
| 6 |
| </p> |
| </td> |
| <td> |
| <p> |
| 1.428456135094165802001953125 |
| </p> |
| </td> |
| <td> |
| <p> |
| 9.41e-007 |
| </p> |
| </td> |
| </tr> |
| <tr> |
| <td> |
| <p> |
| 53 |
| </p> |
| </td> |
| <td> |
| <p> |
| Win32, VC++ 7.1 |
| </p> |
| </td> |
| <td> |
| <p> |
| 13 |
| </p> |
| </td> |
| <td> |
| <p> |
| 6.024680040776729583740234375 |
| </p> |
| </td> |
| <td> |
| <p> |
| 3.23e-016 |
| </p> |
| </td> |
| </tr> |
| <tr> |
| <td> |
| <p> |
| 64 |
| </p> |
| </td> |
| <td> |
| <p> |
| Suse Linux 9 IA64, gcc-3.3.3 |
| </p> |
| </td> |
| <td> |
| <p> |
| 17 |
| </p> |
| </td> |
| <td> |
| <p> |
| 12.2252227365970611572265625 |
| </p> |
| </td> |
| <td> |
| <p> |
| 2.34e-024 |
| </p> |
| </td> |
| </tr> |
| <tr> |
| <td> |
| <p> |
| 116 |
| </p> |
| </td> |
| <td> |
| <p> |
| HP Tru64 Unix 5.1B / Alpha, Compaq C++ V7.1-006 |
| </p> |
| </td> |
| <td> |
| <p> |
| 24 |
| </p> |
| </td> |
| <td> |
| <p> |
| 20.3209821879863739013671875 |
| </p> |
| </td> |
| <td> |
| <p> |
| 4.75e-035 |
| </p> |
| </td> |
| </tr> |
| </tbody> |
| </table></div> |
| </div> |
| <br class="table-break"><p> |
| Finally note that the Lanczos approximation can be written as follows by removing |
| a factor of exp(g) from the denominator, and then dividing all the coefficients |
| by exp(g): |
| </p> |
| <p> |
| <span class="inlinemediaobject"><img src="../../equations/lanczos7.svg"></span> |
| </p> |
| <p> |
| This form is more convenient for calculating lgamma, but for the gamma function |
| the division by <span class="emphasis"><em>e</em></span> turns a possibly exact quality into |
| an inexact value: this reduces accuracy in the common case that the input is |
| exact, and so isn't used for the gamma function. |
| </p> |
| <h5> |
| <a name="math_toolkit.lanczos.h4"></a> |
| <span class="phrase"><a name="math_toolkit.lanczos.references"></a></span><a class="link" href="lanczos.html#math_toolkit.lanczos.references">References</a> |
| </h5> |
| <div class="orderedlist"><ol class="orderedlist" type="1"> |
| <li class="listitem"> |
| <a name="godfrey"></a>Paul Godfrey, <a href="http://my.fit.edu/~gabdo/gamma.txt" target="_top">"A |
| note on the computation of the convergent Lanczos complex Gamma approximation"</a>. |
| </li> |
| <li class="listitem"> |
| <a name="pugh"></a>Glendon Ralph Pugh, <a href="http://bh0.physics.ubc.ca/People/matt/Doc/ThesesOthers/Phd/pugh.pdf" target="_top">"An |
| Analysis of the Lanczos Gamma Approximation"</a>, PhD Thesis November |
| 2004. |
| </li> |
| <li class="listitem"> |
| Viktor T. Toth, <a href="http://www.rskey.org/gamma.htm" target="_top">"Calculators |
| and the Gamma Function"</a>. |
| </li> |
| <li class="listitem"> |
| Mathworld, <a href="http://mathworld.wolfram.com/LanczosApproximation.html" target="_top">The |
| Lanczos Approximation</a>. |
| </li> |
| </ol></div> |
| </div> |
| <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> |
| <td align="left"></td> |
| <td align="right"><div class="copyright-footer">Copyright © 2006-2010, 2012-2014 Nikhar Agrawal, |
| Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert |
| Holin, Bruno Lalande, John Maddock, Johan Råde, Gautam Sewani, Benjamin Sobotta, |
| Thijs van den Berg, Daryle Walker and Xiaogang Zhang<p> |
| Distributed under the Boost Software License, Version 1.0. (See accompanying |
| file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) |
| </p> |
| </div></td> |
| </tr></table> |
| <hr> |
| <div class="spirit-nav"> |
| <a accesskey="p" href="relative_error.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../backgrounders.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="remez.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> |
| </div> |
| </body> |
| </html> |
