A Log Gamma Function Approximation

Daniel D. Warner
June 9, 1999

The Cody-Stoltz Approximation to the Log Gamma Function

The following approximation to the log of the Gamma Function was developed by W. J. Cody and L. Stoltz at the Argonne National Laboratory in 1988. This routine is part of W. J. Cody's specfun package. For background information see the following references.


1. Chebyshev Approximations for the Natural Logarithm of the Gamma Function by W. J. Cody and K. E. Hillstrom, Math. Comp. 21, 1967, pp. 198-203.

2. DGAMMA/DLGAMA, ANL/AMD Program ANLC366S, by K. E. Hillstrom, May, 1969.


3. Computer Approximations by John F. Hart, E. W. Cheney, Charles L. Lawson, Hans J. Maehly, Charles K. Mestenyi, John R. Rice, Henry G. Thacher, Jr., and Chritoph Witzgall. Published by John Wiley & Sons, Inc. 1968.

> Digits := 22;
XINF := 9.99999999999999*10^99; # lambda for FX 2.0
EPS := 8.0*10^(-14); # Machine epsilon for FX 2.0 is not properly defined.
XBIG := fsolve(lnGAMMA(x) = XINF, x);
FRTBIG := sqrt(sqrt(XBIG));
PNT68 := 0.6796875;

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

> lgam1 := proc( xx )
local p, q, d, x, xn, xd, k, z;
# Almost best Uniform Approximation of Degree 7 over Degree 8 on [0.5 - 1.5]
# Note that we assume that xx has already been translated by -1 so that -0.5 <= xx <= 0.5.
p := [ 4.945235359296727046734888, 201.8112620856775083915565, 2290.838373831346393026739,
11319.67205903380828685045, 28557.24635671635335736389, 38484.96228443793359990269,
26377.48787624195437963534, 7225.813979700288197698961 ];
q := [ 67.48212550303777196073036, 1113.332393857199323513008, 7738.757056935398733233834,
27639.87074403340708898585, 54993.10206226157329794414, 61611.22180066002127833352,
36351.27591501940507276287, 8785.536302431013170870835 ];
d := -0.5772156649015328605195174;
x := evalf( xx );
xn := 0;
xd := 1;
for k from 1 to 8 do
xn := xn*x + p[k];
xd := xd*x + q[k];
od;
z := x*(d + x*(xn/xd));
RETURN( evalf(z) );
end:

> lgam2 := proc( xx )
local p, q, d, x, xn, xd, k, z;
# Almost best Uniform Approximation of Degree 7 over Degree 8 on [1.5 - 4]
# Note that we assume that xx has already been translated by -2 so that -0.5 <= xx <= 2.
p := [ 4.974607845568932035012064, 542.4138599891070494101986, 15506.93864978364947665077,
184793.2904445632425417223, 1088204.769468828767498470, 3338152.967987029735917223,
5106661.678927352456275255, 3074109.054850539556250927 ];
q := [ 183.0328399370592604055942, 7765.049321445005871323047, 133190.3827966074194402448,
1136705.821321969608938755, 5267964.117437946917577538, 13467014.54311101692290052,
17827365.30353274213975932, 9533095.591844353613395747 ];
d := 0.4227843350984671393993777;
x := evalf( xx );
xn := 0;
xd := 1;
for k from 1 to 8 do
xn := xn*x + p[k];
xd := xd*x + q[k];
od;
z := x*(d + x*(xn/xd));
RETURN( evalf(z) );
end:

> lgam3 := proc( xx )
local p, q, d, x, xn, xd, k, z;
# Almost best Uniform Approximation of Degree 7 over Degree 8 on [4 - 12]
# Note that we assume that xx has already been translated by -4 so that 0 <= xx <= 8.
p := [ 14745.02166059939948905062, 2426813.369486704502836312, 121475557.4045093227939592,
2663432449.630976949898078, 29403789566.34553899906876, 170266573776.5398868392998,
492612579337.7430887588120, 560625185622.3951465078242 ];
q := [ 2690.530175870899333379843, 639388.5654300092398984238, 41355999.30241388052042842,
1120872109.616147941376570, 14886137286.78813811542398, 101680358627.2438228077304,
341747634550.7377132798597, 446315818741.9713286462081 ];
d := 1.791759469228055000094023;
x := evalf( xx );
xn := 0;
xd := -1;
for k from 1 to 8 do
xn := xn*x + p[k];
xd := xd*x + q[k];
od;
z := d + x*(xn/xd);
RETURN( evalf(z) );
end:

> lgam4 := proc( xx )
local c, x, xs, k, z;
# Assymptotic Expansion with 7 terms on [12 - INF]
# Note that we assume that xx has already been translated by -4 so that 0 <= xx <= 8.
c := [-0.001910444077728, 0.00084171387781295, -0.0005952379913043012,
0.000793650793500350248, -0.002777777777777777777778, 0.083333333333333333333333,
0.0057083835261 ];

x := evalf( xx );
xs := x*x;
z := c[7];
for k from 1 to 6 do
z := z/xs + c[k];
od;
z := z/x;
RETURN( evalf(z) );
end:

> LGamma1 := proc( xx )
local x, z;

x := evalf( xx );
if (x <= 0) or (XBIG <= x) then
RETURN( XINF );
elif (x <= EPS) then
RETURN( -log(x) );
elif (x < 0.5) then
z := -log(x) + lgam1(x);
elif (x < PNT68) then
z := -log(x) + lgam2(x-1);
elif (x < 1.5) then
z := lgam1(x-1);
elif ( x <= 4) then
z := lgam2(x-2);
elif (x < 12) then
z := lgam3(x-4);
elif (x < FRTBIG) then
z := lgam4(x) + (1/2)*log(2*Pi) + (x - 1/2)*log(x) - x;
else # FRTBIG <= x
z := (1/2)*log(2*Pi) + (x - 1/2)*log(x) - x;
fi;
RETURN( evalf( z ) );
end:

The following approximation was included in the Casio Technical Report on the ZX-938 dated April 23, 1997.
This is the Lanczos algorithm that is presented in the book Numerical Recipes in C.

> LGamma2 := proc( xx )
local c, x, t, y, k, z;
# loggamma
c := [ 76.18009172947146, -86.50532032941677, 24.01409824083091,
-1.231739572450155, 0.001208650973866179, -0.000005395239384953 ];

x := evalf( xx );
t := (x+5.5) - (x+0.5)*log(x+5.5);
z := 1.000000000190015;
y := x + 1;
for k from 1 to 6 do
z := z + c[k]/y;
y := y + 1;
od;
# z := (x+0.5)*log(x+5.5) - (x+5.5) + log(2.50662827463100005*z/x);
z := (x+0.5)*log(x+5.5) - (x+5.5) + log(sqrt(2*Pi)*z/x);
RETURN( evalf(z) );
end:

> Digits:=20;lnGAMMA(1);lnGAMMA(1.5);lnGAMMA(2);LGamma2(1);LGamma2(1.5);LGamma2(2);

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

Error Calculation

Since this procedure uses different approximations for several different intervals, it is useful to first look at the absolute error on each interval and then assess the relative error on the entire range. Note that the approximation on the interval from 4 to 12 has the largest relative error. However the maximum relative error is well under the machine precision for the FX 2.0.

> plot((lnGAMMA - LGamma1), 0..0.5, title="Absolute Error");

[Maple Plot]

> plot((lnGAMMA - LGamma1), 0.5..1.5, title="Absolute Error");

[Maple Plot]

> plot((lnGAMMA - LGamma1), 1.5..4, title="Absolute Error");

[Maple Plot]

> plot((lnGAMMA - LGamma1), 4..12, title="Absolute Error");

[Maple Plot]

> plot((lnGAMMA - LGamma1), 12..25, title="Absolute Error");

[Maple Plot]

> plot(((lnGAMMA - LGamma1)/lnGAMMA), 0..25, title="Relative Error");

[Maple Plot]

Now for an examination of the Lanczos approximation, LGamma2(x).
The following results show that the approximation is not quite adequate for the FX 2.0 for small values of x,
and that both the absolute and relative error grow excessively for large values of x.

> plot((lnGAMMA - LGamma2), 0..5, title="Absolute Error");

[Maple Plot]

> plot((lnGAMMA - LGamma2), 0..15, title="Absolute Error");

[Maple Plot]

> plot(((lnGAMMA - LGamma2)/lnGAMMA), 0..25, title="Relative Error");

[Maple Plot]

> plot((lnGAMMA - LGamma2), 0..500, title="Absolute Error");

[Maple Plot]

> plot(((lnGAMMA - LGamma2)/lnGAMMA), 0..500, title="Relative Error");

[Maple Plot]

>