mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2025-11-04 05:50:31 -05:00
578 lines
24 KiB
HTML
578 lines
24 KiB
HTML
|
<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.79.1">
|
||
|
|
<link rel="home" href="../index.html" title="Math Toolkit 2.5.1">
|
||
|
|
<link rel="up" href="../backgrounders.html" title="Chapter 17. 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 17.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 17.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, Jeremy Murphy, 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>
|