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
Following are several tables of coefficients that have been published in different places on the internet.
g=5, n=7 (also called "g=5, n=6")
|
g=7, n=9
|
g=4.7421875, n=15
|
Sample Code
All of the Lanczos tables above are used with the following code (here shown as if written in C):
#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.918938533204672742; 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]; 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); } 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
|
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.
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.
Robert Munafo's home pages on HostMDS © 1996-2012 Robert P. Munafo. about contact
Google+
mrob27
@mrob_27
This work is licensed under a Creative Commons
Attribution-NonCommercial 2.5 License. Details here
s.13