Coefficients for the Lanczos Approximation to the Gamma Function
This page gives several commonly-used sets of values for the coefficients used in the Lanczos approximation of the Gamma function.
Contents
Sample Code
Background
The Gamma function is equivalent to the factorial for integer arguments, but is valid for most real arguments (negative integers and zero being the only exceptions). There are a number of ways to compute it, but each has flaws. Some methods don't converge quickly enough for reasonably small arguments, and some (such as Stirling's approximation actually diverge if you evaluate too many terms. The number of terms that you need to use varies with the value of the argument z and depends on how much precision you need.
The Lanczos approximation has fewer of these problems, provided that the coefficients are computed accurately enough. The Boost documentation14 summarizes the situation thus:
[...] for the Gamma function over real numbers (as opposed to complex ones) the Lanczos approximation [...] appears to offer no clear advantage over more traditional methods such as Stirling's approximation. Pugh12 [...] discovered that [Lanczos and Stirling] were 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:
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.
— John Maddock et al., online documentation for Boost14 library implementation
The Lanczos coefficients are a bit tricky to compute, and the quality of the resulting calculations depends on your choice of two arbitrary parameters g and n. Fortunately, the coefficients only need to be computed once, using a g and n that have been determined (through a thorough search) to be optimal for a given desired precision. The different sets of coefficients given here provide different amounts of precision.
Following are several tables of coefficients that have been published in different places on the internet.
This first set is shown in the sample source code below:
g=5, n=7 (also called "g=5, n=6")
|
g=7, n=9
|
g=9, n=10
|
g=4.7421875, n=15
|
Generating the Coefficients in Maxima
Raymond Toy provided the following Maxima implementation of the formulas by Paul Godfrey3:
load("diag"); Dc(n) := diag(makelist(2*double_factorial(2*k-1),k,0,n)); cmatrix_element[row,col]:= if is(col>row) then 0 elseif row=1 and col=1 then 1/2 else (-1)^(row+col)*4^(col-1)*(row-1)*(row+col-3)!/(row-col)!/(2*col-2)!; C(n) := genmatrix(cmatrix_element, n+1); f(g,n):=sqrt(2)*(%e/(2*(n+g)+1))^(n+1/2); Dr(k) := diag(append([1],makelist(-(2*n+2)!/(2*n!*(n+1)!),n,0,k-1))); bmatrix_element[row,col] := if row = 1 then 1 elseif is(row > col) then 0 else (-1)^(col-row)*binomial(col+row-3,2*row-3); B(k) := genmatrix(bmatrix_element,k+1); lanczos_coeff(g, n) := block([M : (Dr(n) . B(n)) . (C(n) . Dc(n)), f : transpose(matrix(makelist(f(g,k), k, 0, n)))], (M . f)); tlg1(g,n) := float(bfloat(lanczos_coeff(g,n-1)*exp(g)/sqrt(2*%pi)));This sample Maxima session shows the calculation of the coefficients in two of the tables above:
Maxima 5.26.0 http://maxima.sourceforge.net using Lisp SBCL 1.0.57 Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) batchload("toy-lanczos.mac"); (%o1) /home/sample/maxima-work/toy-lanczos.mac (%i2) fpprec:40; (%o2) 100 (%i3) tlg1(5,7); [ 1.000000000190015 ] [ ] [ 76.18009172947146 ] [ ] [ - 86.50532032941676 ] [ ] (%o3) [ 24.01409824083091 ] [ ] [ - 1.231739572450155 ] [ ] [ .001208650973866179 ] [ ] [ - 5.395239384953129e-6 ] (%i4) tlg1(607/128,15); [ .9999999999999971 ] [ ] [ 57.15623566586292 ] [ ] [ - 59.59796035547549 ] [ ] [ 14.13609797474175 ] [ ] [ - .4919138160976202 ] [ ] [ 3.3994649984811887e-5 ] [ ] [ 4.6523628927048577e-5 ] [ ] (%o4) [ - 9.837447530487956e-5 ] [ ] [ 1.5808870322491247e-4 ] [ ] [ - 2.102644417241049e-4 ] [ ] [ 2.1743961811521265e-4 ] [ ] [ - 1.643181065367639e-4 ] [ ] [ 8.441822398385275e-5 ] [ ] [ - 2.619083840158141e-5 ] [ ] [ 3.6899182659531626e-6 ]Sample Code
All of the Lanczos tables above are used with the following code (here shown as if written in C). It is shown with the g=5, n=6 coefficients, but you may instead substitute any of the other sets of coefficients listed above.
// The following constants LG_g and LG_N are the "g" and "n" parameters // for the table of coefficients that follows them; several alternative // coefficients are available at mrob.com/pub/ries/lanczos-gamma.html #define LG_g 5.0 // Lanczos parameter "g" #define LG_N 6 // Range of coefficients i=[0..N] const double lct[LG_N+1] = { 1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5 }; const double ln_sqrt_2_pi = 0.91893853320467274178; const double g_pi = 3.14159265358979323846; // Compute the logarithm of the Gamma function using the Lanczos method. double lanczos_ln_gamma(double z) { double sum; double base; double rv; int i; if (z < 0.5) { // Use Euler's reflection formula: // Gamma(z) = Pi / [Sin[Pi*z] * Gamma[1-z]]; return log(g_pi / sin(g_pi * z)) - lanczos_ln_gamma(1.0 - z); } z = z - 1.0; base = z + LG_g + 0.5; // Base of the Lanczos exponential sum = 0; // We start with the terms that have the smallest coefficients and largest // denominator. for(i=LG_N; i>=1; i--) { sum += lct[i] / (z + ((double) i)); } sum += lct[0]; // This printf is just for debugging printf("ls2p %7g l(b^e) %7g -b %7g l(s) %7g\n", ln_sqrt_2_pi, log(base)*(z+0.5), -base, log(sum)); // Gamma[z] = Sqrt(2*Pi) * sum * base^[z + 0.5] / E^base return ((ln_sqrt_2_pi + log(sum)) - base) + log(base)*(z+0.5); } // Compute the Gamma function, which is e to the power of ln_gamma. double lanczos_gamma(double z) { return(exp(lanczos_ln_gamma(z))); }
Following are lower-precision versions of values in the above tables.
g=7, n=9
|
Issues at Higher Precision
The developers of the Boost libraries describe14 their efforts to implement a numerically useful Lanczos Gamma. They followed the methods of Pugh13 and discovered that it required "an excessive number of terms" to achieve the precision of IEEE 64-bit floating point or higher. For 128-bit precision (such as the IEEE binary128 Quadruple-precision floating-point format) they were unable to find any value of g and n that gave acceptable results.
To handle this problem they transformed the summation at the center of the Lanczos computation (a sum of "partial fractions") into a ratio of two polynomials. This doubles the number of needed coefficients, but also makes all the coefficients have the same sign. This has the great advantage of removing the problem of terms with opposite sign largely cancelling each other out during the summation. They found suitable coefficients with n=13 for IEEE 64-bit, n=17 for 80-bit extended ("long double" on the Intel x86 architecture), and n=24 for a quad-precision HP Alpha format similar to IEEE binary128.
footnotes
1 : Mihai Preda, Implementing the lgamma() function in Java (web article), 2006 Nov 5. This includes an implementation of lgamma from SUN's FDLIBM 5.3, which gets the last few bits of precision for IEEE binary64 (as compared to the 9-term and 15-term Lanczos versions). Java source is here: Preda Gamma.java
2 : Wikipedia, Lanczos approximation. Includes source code of a g=7, n=9 version from the GNU Scientific Library.
3 : Paul Godfrey, A note on the computation of the convergent Lanczos complex Gamma approximation (web page), 2001.
Another version is found on Numericana. One version describes a "15 term expansion" in the abstract; the other a "11 term expansion".
4 : Ken Takusagawa, source code for Riemann-Siegel Z function. Gamma function is in file "cheby.el.cc", function "log_gamma" with coefficients in file "c-coefficients-for-cheby.h", global "const dfloat cof[6]".
5 : Section 6.1 (pp. 213-214) from: William H. Press et al., Numerical recipes in C: the art of scientific computing (2nd edition). Cambridge University Press, (1992) ISBN 0521431085.
6 : Paul Godfrey, Special Functions math library (Matlab code), 2001.
7 : Christian Borgelt, gamma.c (source code), 2008.
8 : Rosetta Code, Gamma function.
9 : Numerical Recipes Forum, Bug in 2nd edition version of gammaln() (discussion thread), July 2005.
10 : Viktor Toth, The Lanczos approximation of the Gamma function (web page), 2005.
11 : Tom Minka, C implementations of useful functions
12 : Glendon Ralph Pugh, An Analysis of the Lanczos gamma approximation (PhD thesis), University of British Columbia, 2004.
13 : Peter Luschny, Approximation Formulas for the Factorial Function (web page), 2011. This page discusses some Gamma formulas, including Lanczos, within the context of approximating the factorial function for large positive real arguments.
14 : John Maddock et al., The Lanczos Approximation as implemented in the Boost C++ Libraries. The C source code can be downloaded, and the implementation of Lanczos Gamma are in files containing "lanczos" and "gamma" in their names, e.g. special_functions/gamma.hpp and special_functions/lanczos.hpp.
Old locations:
https://www.boost.org/doc/libs/release/libs/math/doc/sf_and_dist/html/math_toolkit/backgrounders/lanczos.html
www.boost.org/doc/libs/1_44_0/libs/math/doc/sf_and_dist/html/math_toolkit/backgrounders/lanczos.html
Location that never worked:
https://www.boost.org/doc/libs/1_70_0/libs/math/doc/sf_and_dist/html/math_toolkit/backgrounders/lanczos.html
This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2019 Apr 28.
s.27